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ABSTRACT 
‘ 

A procedure is presented for the analysis of classical subsonic, 
torsion-bending flutter of propellers, extending the earlier work 
by Loring for nonrotating airfoils. Generalized coordinates and 
matrix notation are employed in the analysis; vibration modes 
of a nonrotating blade are employed as basic displacement shapes. 
Making certain simplifying assumptions in order to avoid enter- 
ing into a discussion of the complexities of the most general case 
of coupled engine-propeller system vibration, it is found that the 
flutter equations may be formulated for a single blade without 
loss of generality in considering the mechanism of flutter. 

Because of initial twist, the blade elements undergo vibratory 
displacements in the direction parallel to the relative wind, 
creating additional lift forces and moments, which have been 
incorporated in the analysis. Theoretical expressions are derived 
for the required aerodynamic coefficients. 

An iterative procedure is outlined for the calculation of flexural 
modes of a nonrotating blade. Because of the presence of centrif- 
ugal coupling terms in the equations of equilibrium, changes in 
mode shape with rotational speed are taken into account. 


INTRODUCTION 


er EXTENSIVE AMOUNT of experimental data 
on self-excited vibration of propeller blades has 


been accumulated.'~* Since all of the reported observa- 
tions pertain to blades of wood or improved wood con- 
struction, the problem of flutter prevention has not been 
a matter for serious concern in the design of metal blades. 
However, with the advent of the gas turbine, it is to be 
expected that the need for thinner blades will become 
more acute,‘ and it is not unlikely that the frequency 
spectrum requirements for ftutter prevention may ulti- 
mately establish the limit that can be attained in re- 
ducing blade thickness. 

It is the object of the present paper to describe a 
method of analysis for classical torsion-bending flutter 
of propeller blades. The derivation to be presented is a 
direct extension, employing generalized coordinates and 
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matrix notation, of the earlier work by Loring®® on the 
flutter analysis of nonrotating systems. This treatment 
is believed to be somewhat more general than previous 
formulations of the flutter problem for rotating air- 
foils.2% Attention is confined to the subsonic flutter 
problem, which is likely to be serious, even for a pro- 
peller in which the blades normally operate at supersonic 
speeds. 

The first step in an analysis of this kind is to replace 
the continuous elastic body by a mathematical idealiza- 
tion with a finite number of degrees of freedom. This 
objective is accomplished by the selection of a set of 
basic displacement forms. The allowable displacements 
of the idealized system then consist of linear combina- 
tions of the basic shapes, with coefficients that are 
functions of time (generalized coordinates). Lagrange’s 
equations are employed to establish the equations of 
motion relative to a rotating coordinate system, with 
generalized inertia, elastic, centrifugal, and aerodynamic 
forces included. 

In the flutter analysis of wings that carry no large 
concentrated masses, the principal requirement is that 
the basic shapes shall be sufficient to represent (in suit- 
able linear combinations) the first two or three normal 
modes of the wing. It is proposed here to employ 
vibration modes of the nonrotating propeller as basic 
shapes. The influence of the centrifugal field of force 
on vibration mode shapes is approximately taken into 
account (if a sufficient number of modes of the non- 
rotating blade are included in the analysis) by the pres- 
ence of certain centrifugal coupling terms in the equa-' 
tions of motion. 

The following theory is derived for blades whose mass 
centroids lie on the elastic axis; a trivial extension 
would be required to treat the more general case. 
Since the amount of initial twist is small in a normal 
blade, it is concluded that the elastic coupling between 
bending and twisting may be neglected. The elastic 
axis is assumed to be initially straight and perpendic- 
Simplification of the 


ular to the axis of rotation. 
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general problem in order to confine consideration of the 
mechanical system to the propeller and to avoid the 
mathematical complexities of coupled engine-propeller 
vibration is discussed in Appendix B, which proposes 
limiting the analysis to a propeller whose hub is con- 
strained against lateral and axial vibration and is not 
subject to torque reactions from an attached power- 
plant or drive system. 

A preliminary vibration analysis of the nonrotating 
blade is carried out to determine the first torsional 
mode and corresponding natural frequency. Ordinarily, 
it will be desirable to include in the analysis all of the 
bending modes with natural frequencies less than the 
first torsional frequency. With regard to boundary 
conditions, under the assumption of the isolated pro- 
peller already mentioned, it will be necessary to con- 
sider both hinged and clamped support conditions for 
the blades (resonant and antiresonant modes, respec- 
tively). Mode shapes for the entire propeller can be 
obtained by combining mode shapes of the individual 
blades in such a way that the boundary conditions are 
satisfied for the propeller as a whole. 

It is shown in Appendix B that resonant and anti- 
resonant modes are uncoupled and that the flutter 
equations may be formed for a single blade. Until the 
time when sufficient experimental work with modern, 
high-performance propellers will have indicated the 
types of modes to be expected in potential propeller 
flutter problems, hence the nature of boundary condi- 
tions leading to the most critical flutter conditions, it 
seems advisable to make separate investigations into the 
possibilities of flutter in resonant modes and in anti- 
resonant modes. An analogous situation occurs in the 
case of wing flutter, where both symmetrical and anti- 
symmetrical flutter are usually investigated. 

Perhaps the greatest obstacle to the analytical treat- 
ment of propeller flutter is the lack of a sufficiently 
general airfoil theory. In the following developments, 
theoretical lift and moment coefficients for an oscillat- 
ing airfoil of infinite aspect ratio in a uniform rectilinear 
flow are used as section coefficients. Some of the ob- 
vious objections to this procedure are as follows: (1) 
Three-dimensional considerations, such as tip effects, 
helicoidal form of trailing vortex sheets, etc., may be 
significant; (2) blade interference effects may be im- 
portant; (3) oscillatory drag forces may become im- 
portant when the blades are operating at high angles 
of attack. These objections appear least serious for 
propellers with blades of high aspect ratio operating at 
low angles of attack and high values of the advance 
ratio. Because of initial twist, the elements of a 
vibrating blade exhibit significant components of dis- 
placement in the direction of the relative wind, and 
additional aerodynamic forces are created. Contribu- 
tions to lift and moment due to motions of this type are 
included in the present scheme of analysis. Formulas 
for the required aerodynamic coefficients are derived in 
Appendix A. 
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FORMATION OF THE FLUTTER EQUATIONS 


Lagrange’s equations for a dynamical system with a 
finite number of degrees of freedom are to be used in 
forming the flutter equations. For reasons presented 
in Appendix B it is permissible to deal with a single 
blade in forming these equations. 

For motions of small amplitude the apparent kinetic 
energy of the vibrating blade, relative to the rotating 
reference frame, may be represented by a quadratic 
form in the generalized velocities 


T = ('/2)I,q’aq (1 
where 
g = column of generalized velocities 
a@ = nondimensional inertia matrix 
I, = reference parameter with dimensions of a mo- 


ment of inertia, lb.-sec.?-in. 


+ 
The potential energy is represented by a quadratic 
form in the generalized coordinates 


V = ('/2)K,g'(f + fea \¢) 


~ 


where 
q = column of generalized coordinates 
f = nondimensional elastic stiffness matrix 
fe = nondimensional centrifugal stiffness matrix 


K, = reference stiffness, in.lb. 


The representation [Eq. (2)] encompasses two forms of 
potential energy: that arising from the work done in 
opposition to internal elastic forces and that done in 
opposition to the centrifugal field. 

The column of generalized aerodynamic forces, due to 
steady sinusoidal oscillation, may be expressed in the 
form 

Q = —pwby*Rcq (3) 
where 

p = air density, lb.-sec.*/in.* 

: angular frequency of oscillation, rad. per sec. 
bo = reference semichord length, in. 

R = tip radius, in. 

¢ = nondimensional air-force matrix 

For steady oscillations the equations of motion are of 


form 
[—I,w*a + w*mpbo'Re + K,(f +fedl¢g = 0 (4) 


Explicit formulas for the matrices a, c, f, and f, will be 
derived in the following sections. It will be found that 
the elements of f, are proportional to the square of the 
rotation speed. Furthermore, the elements of ¢ are 
functions of reduced frequency, Mach Number, and 
advance ratio. 

For computational convenience Eq. (4) is multiplied 
by 


p?/w?arpbotR 
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where 
u? = ko*/(1 + ho’) 
ko = reduced frequency at reference station = who/ Vo 
lV) = resultant air speed relative to blade at reference 


station 
Then it follows that 
lo + X(f + fd]Jg = 0 (5) 

where 

( ‘u(c — Ta) 

i = i mpby'R 

X = p*(a,/w)? = (V,/Vo)2(1 — pw?) 

w, = reference frequency = +/K,/mpbo'R 

V, = reference air speed = w,bo 
The determinantal equation is of form 


Co 4 X(f + Te) = 0 (6) 


METHOD OF SOLUTION 


The stability of the propeller is to be examined at 
several different rotation speeds, 2, while the advance 
ratio and pitch setting are held constant. At each value 
of Q, 
of Eq. (6) are determined. 
of X are obtained, which may be expressed in the 


several values are assigned to uw, and roots 
In general, complex values 


form 
X = X,(1 — 7g) (7) 


The solution of Eq. (6) consists of several distinct, con- 
tinuous functions of uw (branches), each generated by 
one of the roots of this equation. 

In the analysis of nonrotating systems, the generalized 
coordinates are usually chosen so that there is no elastic 
coupling; the diagonal elements of f are of form 
(1 + 2g;)fy, where g; = structural damping factor as- 
sociated with the jth coordinate When gis small, it isa 
measure of the stability of the system, in the sense that 
if all of the structural damping factors were reduced by 
the same amount, equal to g, then a condition of neutral 
stability would exist at the given value of u (reference 
6, pp. 121-122). It is not possible to make precisely 
the same physical interpretation of g for a propeller, 
since the elements of the elastic stiffness matrix occur 
in combination with the elements of f,. A similar 
analytical procedure may be employed in propeller 
analysis, however, and g may be regarded as a measure 
of stability in a slightly generalized sense. 

The frequency and air speed at the reference station, 
corresponding to a particular value of X,, are calculated 
from the formulas 


wie, = w/VX 


V/V, = Vi — w/VXS (9) 


However, since 2 and U//QR are given in advance, the 
following relation must be satisfied: 
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, (w,b9/ QR)? 
Xi = se (1 


— 9 
(U/QR)? + kr 7 9) 


where {9 = 20/Rand z = radius at reference station. 

The correct values of X, and yu at the given value of 
Q are determined by plotting Y; vs. u for each branch 
and locating the intersections with a plot of Eq. (9). 
Auxiliary plots of g vs. u may be used to determine a 
value of g for each branch. By repeating this process 
for additional values of Q, sufficient data will be ob- 
tained to construct plots of g vs. 2 at the assigned value 
of U/QR. Critical speeds occur when g becomes equal 
to zero. 

Because of the fact that the composite stiffness 
matrix (f + f.) is not of diagonal form, there is some 
increase in the computational labor required to solve 
Eq. (6). However, if the number of coordinates is not 
greater than five, the expansion of the determinantal 
equations can be accomplished without much difficulty 
on any good desk-type calculating machine. The La- 
place expansion provides a suitable computational 
scheme for this purpose. 


VIBRATION MODES OF NONROTATING BLADE 


The fundamental torsional mode is required. Several 
well-known methods are available for its determination; 
hence, this phase of the problem will not be discussed 
in the present paper. The theory of flexural vibrations 
for a slightly twisted beam is discussed in some detail, 
and a procedure is presented for the calculation of anti- 
resonant modes. A similar scheme may be employed 
in the determination of resonant modes, with suitable 
modifications to account for the absence of hub re- 
straint. 

The orthogonal axes—OX, in the plane of rotation, 
OY, the axis of rotation, and OZ, a radial axis along 
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the blade—are introduced with the basic assumption 
that the center of mass and the elastic center are co- 
incident and that the locus of these points is a straight 
line coinciding with the Z-axis. The intersection of the 
Z-axis with a plane parallel to the X-}-plane estab- 
lishes the local origin, O’, for the £»-coordinates, corre- 
sponding to the minor and major axes of inertia of the 
blade section. £8 is the local angle of twist measured 
from the X-Z-plane. 

The typical blade element is subjected to a system of 
forces and moments parallel to the X-Z-plane, as shown 
in Fig. 2. Displacement, shear, and bending moment 
in this plane are identified by the subscript 7, suggested 
by the radial character of the centrifugal field in this 
plane. 

The equations of equilibrium are 


> Fz: dT/dz = —Q?mz (10) 
> Fy: dS,/dz = — (w? + 2*)mu, — (d/dz)(Tdu,/dz) 
(11) 
> M: dM,/dz = —S, (12) 
where 
m = mass per unit length of blade, lb.-sec.?/in.? 
T = centrifugal tension, Ibs. 


S, = shear force, lbs. 


M, = bending moment, in.Ibs. 

# = angular frequency of vibration, rad. per sec. 
Q = rotational frequency, rad. per sec. 

u, = vibratory displacement in direction of X-axis 
6 = inclination of elastic axis = du,/dz 


A similar analysis for the ¥-Z-plane leads to the equa- 


tions 
> Fy: dS,/dz = 
>: dM,/dz = —S, 


(d/dz)(Tdu,/dz) —_ (13) 


—w"MUy, 
(14) 


It is noted in the derivation of these equations that the 
centrifugal force has no component in the direction of 
the Y-axis. The identifying subscript p is suggested 
by the parallel character of the centrifugal field. 
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By integration of Eq. (10) there is obtained for the 


centrifugal tension 
R 
f mQ*2dz 15) 


From Eqs. (11) and (12) it follows that (primes de- 
noting differentiation with respect to z) 


T = 


M,” = m(w? + Q?)u, + (Tu,’)’ (16 
similarly from Eqs. (13) and (14), 
M,” = 


mat, + (Tuy’)’ ~ (17) 


By integration of Eqs. (16) and (17) and introduction 
of the boundary conditions 1/, = MV,’ = M, = M,' = 
T = Oatz = R, it follows that 


-R PR oR 

M, (w? + x) | fm dadz — f twa: IS 
2 PR oR 

M, = of J mupteds - | Tu,‘dz 19 


If the notation is condensed, according to the scheme 
introduced by Arnoldi,® by introduction of the column 





matrices 
{Mp} = {M,, M,}, {trp} = {tr Up} 
diagonal matrices 
[m(sz) O 
[m] = 
LO m(z) 
Pw? +2? 0 
[u?] = 
| 0 w* 
T(z) O 
[7] = 
GO zits) 
and matrix operators 
¢ O 6 O 
[o] = » [6] = 
0 @ 0 6 
where 
*R 
f(s) = df/dz, of(z) = f(z)dz 


then Eqs. (18) and (19) may be reduced to a single matrix 
equation 


f eee l 
{Mu} = 


[w?] [o]?[m] {rp} — [o}[T] [6] {up} (20 


Simple bending theory provides the relation .1/ 


EIx for an untwisted beam, where 


FE = modulus of elasticity, lbs. per sq.in. 
I = principal moment of inertia, in.* 
k = curvature in the principal plane, in.~' 


In applying the theory to a twisted propeller blade, it 
will be necessary to introduce the transformation from 
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the £-n-coordinates to the x-y-system. With reference 


to Fig. 1, 
j Bea f \ 
1 Men j —_ [8B], Urp} 
where 
fun} = {uz, uy} 
u displacement in direction of £-axis 
u displacement in direction of 7-axis 
cos B sin 8 
3) 


— sins cos~p 


Furthermore, it is easily shown that 


f i f \ 9 
1 Meny = [8B], Mp} (21) 
f \ f l 929 
Ken = [B] rps (22) 
with 
Vi bending moment in the &-s-plane 
kK; curvature in the &-s-plane 
\/, = bending moment in the 7-z-plane 
k, curvature in the -s-plane 


Positive values of xk, .\/; occur where the £-coordinate 
of the center of curvature is positive; similarly, «,, WZ, 
are positive where the y-coordinate is positive. In the 
é-n-coordinates we have the relation 


{Men} = 1B) {xe (23) 
where 
[B] = 
0 EI, 


and from Eqs. (21), (22), and (23) it follows that 
{kp} = 18] -[B]-"[8] { Mop} 
But [k,,] in untwisted coordinates is approximately 
equal to the second derivative of the deflection 
{Kr} = [5]? { 2} 
hence, 
[5]?{u,>} & [8] [B]-"[8] {AL} (24) 


Integrating curvatures to obtain deflections and intro- 
ducing the boundary conditions u, = u,’ = u, = u,’ = 0 
ats = 0, we define the matrix operator 


o 0 
lo] = 
0 o 
where 
af(z) = f(z)dz 
0 


é is an integral operator similar to the o form already 
defined but differs in that curvatures must be integrated 
irom the origin outward to the current station, whereas 
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integrations for moments are performed from the cur- 
rent station to the tip. Finally, the desired relation 
between displacements and moments is 
{ur} = [o]?[8]-'[B]-'[8] { AL} (25) 
Substitution of Eq. (25) into Eq. (20) yields 
Mint = ([w*] [o]?[m][o]* — [o][T][o])[8]-' Xx 
[B]—[8] {Mp} (26) 


an integral equation in operational form which can be 
solved by a variety of methods. An alternative form 
referred to the &n-coordinates may be obtained by 
introducing the transformation [Eq. (21)]; the result is 
{Men} = [8]([o*] [o]?[m] lo]? — [o][T][o]) X 
(8]-"[B]-"} M,,} (27) 
It is convenient at this point to derive certain rela- 
tions that will be needed later in the formation of stiff- 
ness and inertia matrices. The total potential energy 
for the deformed blade is given by 


+ ig y » 
VY = sf M «dz + 5 / M,x,dz = sf Mu,” dz + 


l 
2). 


Integrating by parts twice and introducing the appro- 
priate boundary conditions, we obtain the alternative 


Fe 1 f* 
V = sf M,"u,dz + 5 Motus (28) 
“/0 -J0 


Combining Eqs. (16), (17), and (28) gives 


° a o2 *R 
5 w* * - AY 4 ~ 
V=-— m(u,? + u,*)ds + mu,2dz + 
9 B 9 
0 ~ J 0 


R 
” ” 
M,"u,"dz 


expression 


F i- 
| [(Tu,')'u, + (Tu,')'u,|dz 


which through integration by parts may be converted 
into the more convenient form 


2 R 2 R 

, w 9 9 9 

V = = m(U,? + Uy*)dz + . mu,2dz — 
2/0 Z 0 


Pm 

l kh / 4\9o > 

sf T [(u,’)? + (up’)?]dz (29) 
0 


~*~ 
~ 


The quantity, V, given by Eq. (29) is in reality equal 
to the elastic strain energy that is stored in the deformed 
blade. The sum of the second and third integrals 
(generally negative) represents the work done on the 
blade by the centrifugal field during the process of 
deformation. In forming the flutter equations we shall 
require an expression for the potential energy arising 
from work done in opposition to the centrifugal field; 
the required expression is 


a f® Lf 
i. mu,*dz + 5 T [(u,')*? + (up’)?]dz (30) 
Z/0 -J/0 
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FIG. 3. 


The total potential energy that appears in Lagrange’s 
equations is the sum of V, and the elastic strain 
energy. 

For computational purposes the continuous blade is 
replaced by a set of NV interconnected massless links of 
equal length carrying concentrated masses at the 
hinge points, as indicated in Fig. 3. 

Instead of continuous components of curvature, 
k,(Z), Kp(z), we have a discrete set of angular displace- 
ments between successive segments 


Also, introducing the column matrix 





{Mo} = {M,(20), M,(a), ..., M,(ey 
and the following square matrices (of order 2.V) 
[ #7¢(20) 
E:1,(21) 
EI, (z 
[B] : 
[ ¢08 By 
COS Bo 
COS Bx 
(3) 
—sin ~; 
—sin Be 
len —sin By 





1) 


Azk; : 
Pon ra \ = components of angular displacement of ith 
aii segment relative to (i — 1) st 
The integral operators o and o are redefined as follows: 
N 


of(2) = DU f(z) dz 


© let 


i 
d Sf(z,) Az 
j=l 


af(z) 


Also 
r $..,3 
( scart 
bw eee | 
0 O I 
lo] = [o]’ 
my 
Me 
[m] = 
Mn 
ae ee o..% 
og Vol “2) oN4j 


Then for the segmented blade, we have, instead of 
Eq. (15), 

{T} = Q{o][m]{z} = a{T} 31) 
where 


{T} = [o][m]}z} 32 





1); M,(20), M,(2:)..., M,(ew-1)} 
ETI, (20) 
EI, (21) 
I 4 oN 1) 
sin By 7) 
sin Be 
sin By 
cos §; 
COS Be 
cos sel 
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[co] O(N X N) 
[i] = 
O(N X N) [o] 
[m] O(N X N) 


(S] = (XI 


[M] = 
O(N X N) [m7] 
we obtain as the analog of Eq. (26) for the nonrotating 
blade (JT = 0) 
(Mp} = o[DO)M1S) 218) 1B) 16] { ap} 

or 

[C — rx] { M,,} = 0 (33) 
where 

A = 1/w’? 

[(C] = [A][D]  _ 

[4] = (14) (2)? 

[D] = [8]~"[B]~“[8] 


Iteration procedures may be used to obtain the latent 
roots and modal columns of the matrix [C]; punched 


card machines are particularly suitable for this phase of 
the analysis. This is a familiar problem, and a dis- 
cussion of computational techniques would be out of 
place in the present paper. , However, since [1], [B] 
are of diagonal form and 


[2] = ()’, [8]-* = [8]’ 


it follows that both [A] and [D] are symmetrical. 
Consequently, if {1/,,“} is a modal column corre- 
sponding to the latent root \,, then the modal row is 
equal to 


{Mg }'1D] = {9} 


This property reduces appreciably the labor required 
for the determination of higher modes, since the ortho- 
gonality relations involve both modal rows and modal 
columns. If the blade were treated as a continuous 
system and conventional methods of numerical inte- 
gration (based upon polynomial approximations) were 
employed in forming the matrices [=], [2], then the 
relation [>] = [Z]’ would not hold. It is for this reason 
that the segmented blade has been introduced. 


CENTRIFUGAL STIFFNESS MATRIX 


When the continuous blade is detormed in pure bending, the potential energy due to the opposing action of the 
centrifugal field is given by Eq. (30). The corresponding expression for the segmented blade is 


Ve = (PAS) (_ 4} foe fae} +, { typ’ } [7] { u,’} ) (34) 


» = 


where 


{u,} = { u,(z1), u,(Ze). . .U, (Sy) } = [o]*{x,} 


{ thy’ } = {14,'(20), Uy’ (21),. « « 52y’(Zy—1)3 Up’ (Zo), My’ (Z1) .. 5 Mp (2y—-1)$ = [2] { Kp} 


PT 


T; 


[T] = 





Ty 


T, 


to 





TJ 


For the segmented blade, u,’(z;_ 1), u»’(2;~-1) are the slopes of the 7th segment. The diagonal elements of the matrix 
5 r t—1 7 #—1 g g 


[7] are obtained from Eq. (32). 


The column of generalized coordinates for the flutter 
analysis is denoted by q = { aq, “er In} « The first 
n—1 q's correspond to bending modes; 4, is associated 
with the torsional mode. Displacements and slopes are 
related to the generalized coordinates by matrix equa- 
tions 


{ttt = [Urslg (35a) 
ttrp’| = [U»']@ (35b) 


[U,,] and [U,,’] are rectangular matrices with 2N 
rows and n columns; typical columns of these matrices 
characterizing the 7th bending mode are, respectively, 
f ay . y72h5, © 26. 
Up 5 = [=] { Krp } (36a) 


.,n—1) (36b) 


f Ge a ae. Cee ee a. ‘ 
} Ury’ } -_ [2] } rp ‘y (7. = i. - 4 ° 


The last column of both Eqs. (35a) and (35b) is com- 
posed entirely of zeros. 
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Inserting Eqs. (35a) and (35b) into Eq. (34), we ob- 


tain 
Fe = (} 2) Kig'f eng (37) 
where 
0? Az --w - t r ween , ) 
eer (—[U;] [m][U-] + [Un'IT NU») (28) 


The final row and the final column of Eq. (358) contain 
only zeros. 

A simple development leads to the following expres 
sion for the contribution of the torsional mode to the 


centrifugal energy, Vc: 


09," . 
Va = | a(/; — I,)¢* cos 2pd) (39) 


where 

l, major moment of inertia of blade section = 
S S €dédn, in.! 

i, = minor moment of inertia of blade section = 
S S wdédn, m.* 

Dp = density of blade, Ib.-sec.*,in.* 

¢(z) = function of z representing the first torsional 
mode 


Finally, the matrix f, is obtained by adjoining the di- 


agonal element 
oR 
an K,) f pl. — I,)¢? cos 28dz (40) 
0 


to f, 


INERTIA MATRIX AND ELASTIC STIFFNESS MATRIX 

From the orthogonality of the basic shapes it follows 
that the inertia matrix, a@, and the elastic stiffness 
matrix, f, are of diagonal form. The elements of the 


principal diagonal of a are given by 


dy = (Az/T,) ftp} [M1] fury? } (é = 1,2, .., 2-1) (41) 


and 
"R 
Gin = (1 1) f [,¢°dz (42) 
. 0 
where /, = polar moment of inertia per unit length of 
blade, referred to elastic axis, lb.-sec.2 = p(/; + J,). 


The typical diagonal element of f is given by 


fu = (w?/wo?)ax (4 = 1, 2, ..., 2) (43) 


where 


w, = natural frequency associated with 7th coordi- 
nate at 2 = 0, rad. per sec. 
wo = WK,/I,, rad. per sec. 


AIR-FoRCE MATRIX 


As indicated in the Introduction, aerodynamic coefti- 
cients for two-dimensional flow are to be employed as 
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section coefficients in the calculation of generalized 
aerodynamic forces. With reference to Fig. 4, consider 
an airfoil oscillating harmonically with small amplitude 
and having an initial angle of attack equal to a. The 
displacements to be considered are: 


bl = displacement of one-quarter-chord point in 
opposition to relative wind 


bA = downward displacement of one-quarter-chord 
point in direction perpendicular to relative 
wind 

B = rotation of airfoil, positive in the stalling direc- 
tion 


In the notation of Kiissner and Schwarz™ (with addi- 
tional coefficients k,, m, and 6 in place of /) the lift 
and moment are given by 


K = mrpV°*b[ITk, + Akg + Bho] | 


Mo = mp V°b?[ITm, + Am, + Bm] f (4) 


The available evidence'!'* indicates that k,, &,, mm, 
m, do not vary significantly in the range |a| < 5°. 
Explicit formulas for these coefficients are given by 
Kiissner and Schwarz’? for incompressible flow. Theo- 
retical values incorporating the effects of compressibil- 
ity are available from other sources" for a considerable 
range of Mach Numbers. It is shown in Appendix A 
that 


kn = alka + 2ik//1 — M*)\ 


My, = QM, f 
For incompressible flow it follows that 


kn = al—{k? + 2kG) + 2ik(1 + F)]) 

mM, = —ak?/2 { 
where F(k) and G(k) are the functions introduced by 
Theodorsen.'* The lift acting on an airfoil at fixed 
angle of attack in a stream with harmonically varying 
velocity has been determined by Isaacs.'® Greenberg" 
has also developed formulas for the lift and moment on 
an oscillating airfoil in a pulsating stream. Both of 
these studies are concerned with incompressible flow, 
in which the amplitude of velocity variation need not 
be small in relation to the mean velocity. Eqs. (46) 
are in agreement with Greenberg’s results for the case 
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in which the variation in main stream velocity is of the 
same order of magnitude as the velocities of airfoil 
oscillation. 

From this point the derivation of the airforce matrix 
follows Loring’s development® closely. However, cer- 
tain modifications in detail are required as a conse- 
quence of the variations in air speed and reduced 
frequency along the blade. In the calculation of 
virtual work it is convenient to employ the “‘displace- 


ments” 


hy 


ll 


(b/bo) H 


hy = (b, by) A, he = 5 
and corresponding “‘forces”’ 
M, =0, M. = —dK, M3 = —My 


(Information with regard to the oscillatory drag force 
is not available at present; 1/; is introduced here in 
order to simplify the subsequent derivations.) It will 
be observed that /,, iy represent linear displacements 
referred to bo, instead of the current semichord, b. 
The column of ‘‘forces’’ M = { My, M2, Ms} is related 
to the column of “‘displacements” = {hu, ho, hs} by a 


matrix equation 
M = —pwhpOAGh (47) 


where 


j= b, bo 


0 0 0 
A = ky k? k, k? k,/k? 


9 


m,/k? m,/k? m,/k? 
a 


Furthermore, the column h is related tc the column, q, 
of generalized coordinates by an equation 


h = Lq 


The elements of the matrix L are functions of the radial 
coordinate, 2, representing the basic displacement 
shapes referred to local wind axes. A straightforward 
application of the principal of virtual work leads to 
the following formula for the column, Q, of generalized 


aerodynamic forces: 


Q = —mpwbo'Rcq (48) 


at 
c= f L/O@AOLdté (49) 


where 


f = 2/R 


Numerical difficulties will be encountered if the lower 
limit of integration, (1, is too small. Air forces on the 
blade near the hub may be neglected. 

The resultant air speed at radius g is 


V = VU? + 2%? 


where U = forward speed plus inflow velocity, assumed 
constant, and at the reference station 


Vo = VU? + 2%? 
Hence, the reduced frequency is given by 
k = wh/VU? + 2s? = kb 
where 
6 = (b/bo) V (U/QR)? + §0?/ V (U/QR)? + & (50) 


For incompressible flow the matrix A may be ex- 
pressed in the following form introduced by Loring® 
(Az: may be deleted in the present ’case) 


QF 2G f/1 
) 2F , 
(je a» + ( k )s ie 


where 

rT 0 0 O hy 0 O 

l 
a. = oo 0 1 

ica a l 5 | —_ 0 

a l 3 
— =. 0 0 0] 
0 0 0 0 0 0] 
Ay=la l 1 |, A; =|2a 0 1 
0 0 0 0 0 1 | 











When solutions are required for many different values 
of reduced frequency, a considerable saving in time can 
be accomplished by introducing Loring’s approximate 
representations for F and G as follows: 


J F(k) F (yo) = Fy + ¥1(7 ) Fo’ a Yo(y) Fo” 
\G(k) = G(yko) = Go + valy)Go’ + Yaly)Go" 


Il 


where 


| Fy = F (ko), fF,’ = F(9.40Ro) — fy 
“” = F(2.00ko) — Fo 

Go = G(ko), Go’ = G(0.40k9) — Go 
be: = G(2.00ko) — Go 


(52) 


j@ for a nonrotating airfoil 
(6 for a propeller blade 


The functions y, have been tabulated.® Substitution 
of Eqs. (51) and (52) into Eq. (49) leads to the formula 








& 
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c= [C\ + (2F ko") C3 ot (2F)’ ho?) Cs - (2F 5” ko”)Cs = (2Go 'Ro) Ce —_— (2Go’ ko) C;’ — (2G," ko) 3'| + 
1[(1 ko) Cy + (2Go/ko?) C3 + (2Go’/ko?) Cy’ + (2Go" /Ro?) Cs’ + (2Fo/Ro) Ces + (2F)’ ko) Cr + (2 0” / Ro) Cs] (53) 


where 


a = S Md, 
Cs = J’ra(1/6) ede, 
Cs = SM O)Yde, 
= Sr3(1/8) vsdb, 
Cy = fra(1/6)esde, Cs’ = 


C3 = J’rs(1/6)*de, 
Cs = Sr4(1/0)de, 
Cy = Srs5(1/0)de 
Cs’ = Sr3(1/6) Pade 
S (1/0 


C4 - S r3(1 6) "Wide 
C; = Sl O)yidt 


Ai = L'OAOL, Yi = ¥i(8) 

Finally, the matrix co), which appears in Eq. (6) is given by 

_ 2h Ci. 2F,’ .. 2 Fy” 2Goko 2Go'ko 2G)", 
a = Cy + ~ C3 es =e r . eT ee = Ue se - C7’ = = C’ + 

] + ko l 2 ky? l + Ry? l T ky? l + ko? l + ko? t+ ky? al 
oo. 2Go 2G,’ 2G" 2 Foko 2 Fo’ ko 2F "ko | 
1 Cor ~C3 + -C,' + - C,’ + - aor a C (54 
| + ho’ 1+ ho.) 1+ hy’ lem” ita” 140" ttn 

with C\’ = C, — Ta. The functions of kp appearing in Eq. (54) have also been tabulated.*® 


It is advisable to begin the flutter analysis of a pro- 
peller by surveying the complete range of reduced fre- 
quency for incompressible flow. Subsequent calcu- 
including compressibility effects. may be 
The re- 


lations, 
limited to only a few significant values of po. 
quired air-force matrices may be calculated by direct 
integration for each value of o, in accordance with Eq. 
(49). The radial variation of the Mach Number is 


given by 


/ ¥ : = i Seg FE 
M/Mo V (U/QR)? + §°/V (U/QR)? + $0? 
where 
Jo = Mach Number at reference station (QR/c) X 
V (U/QR)? + Fo? 
( speed of sound in undisturbed stream 


However, in view of the limited data on compressi- 
bility effects which are available at the present time, it 
appears reasonable to treat the Mach Number as a 
constant, equal to the correct value at some typical 
station, say ¢ 0.80. 
APPENDIX A.—LiIFT AND MOMENT DUE TO 
TRANSLATORY OSCILLATION 

We consider an airfoil of infinite aspect ratio at a 
fixed angle of attack, a, executing translatory oscilla- 
tions of small amplitude parallel to the undisturbed 
flow. The airfoil displacement, positive in the direction 


of the negative x-axis, is 
h(t) = bHe™ (Al) 


Relative to fixed x-y-axes, the components of fluid veloc- 
ity are U+ u,v. Referred to &-n-axes moving with the 
airfoil, the velocity components are U 
Any physical quantity associated with the moving 
particles of fluid may be represented by a function 


+h, + u, v. 


f = f(é& 0, 2) 


The same physical quantity referred to fixed axes is 
represented by a function 


F(x, y, t) = f(x + h, y, t) 
It follows that 
4 = tes Fy _ | = A —_ hi fe + Si (A2) 


Relative to the x-y-axes the equations of motion of a 
fluid particle are of form (body forces absent) 


u, t+ (U + u)ur + vu, = —(1/p)p,) 


AS 
y —(1/p) py j 


y.+ (U + uo, + 
However, if u, v, p are given initially as functions of 
£, n, t, then it follows from Eqs. (A2), (A3) that 

ut+(U+h,+ uu; + vu, = —(1/p)p, | 


; Ad 
uy+(U+h,+ u)v, +, = —(1/p)p, § 


4 


By introducing the usual linearizing assumption (4%, 
v < U), Eqs. (A4) are reduced to the form 


— | l p)pe 


n+ (U + iden = —(1/p)p.} 


It is assumed further that the flow is irrotational, 


permitting the introduction of a velocity potential 


(£, 7, ) such that 


u = dz, v= @q, 
Then Eqs. (A5) may be written 
[Io + (U+ hy) de] = —(1/p)p; t \6 
[Io + (U+ hide], = —(1/p)p, f 


By introducing the assumption of adiabatic flow (so 

that the density is a function of pressure only) with no 

disturbances at § = —o and neglecting changes in 
6 


density, Eqs. (A6) may be integrated to obtain the 
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relation 
p/p = C—G& — (UF lid: (A7) 
with C = constant. 
The continuity equation may be written in the form 
uz + v, + (1/p)(dp/dt) = 0 (AS) 


Furthermore 


dp/dt = (dp/dp)(dp/dt) = (1/c*)(dp/dt)  (A9) 


where c = speed of sound in undisturbed flow. From 
Eq. (A7) it follows that 
1/po)(dp/dt) = —du — hide — (U + Inde: — 
(U + hi + u) [be + (U + hide] — 


(A10) 


V[dy + (U + hidden] 


' By inserting Eqs. (A9) and (A10) into Eq. (AS), de- 


leting terms that contain u/c, v/c, h,/c as factors, and 
replacing u, v by ¢;, ¢,, the following differential equa- 
tion is obtained: 

1— W°")dbee + 6, — (2M /c) be — (1/c7)by = 0 (A11) 
where MJ = U/c. 
form as its counterpart in the linearized theory referred 


Since this equation has the same 


to stationary axes, it follows that the solution of the 
present problem may be expressed in terms of results 
that are already known. 

The component of fluid velocity normal to the airfoil 
must vanish. Consequently, 


v= —al(U+h,+ u) (A12 


The term au is neglected, and the boundary condition 
is imposed at the £-axis. Hence, the condition to be 


satisfied is 


—a(U + h), ~6 S$ = & 6, uy] = @ 


Also, it is required that u remain finite at the trailing 
edge E b). 

3ecause of the linearity of the problem it can be 
solved in two parts. We consider first the quasi-steady 
problem with 


v=, = —aU, —bs8 tit), 1 =0 


In this case there is a solution of (A1l1) independent of 
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t. From Eq. (A7) with ¢, = 0, it follows that the pres- 
sure jump, lift, and moment are obtained by introducing 
the factor [1 + (4,/U)] into the formulas for the steady- 
state case. The oscillatory lift and moment are 
K® = Q¢a/W1 — M?)(pU?/2)(2b)(iy/U) = 

rp U%b(2iak/V1 — M?)H (A13) 
M;” = 0 


The second part of the solution is determined by the 
boundary condition 


Oa] 


= v= — ah, 


Hence, 
the term in Eq. (A7) containing /, is of second order; 


¢@ and its derivatives are proportional to /,. 


it will be neglected. The problem then becomes equiv- 
alent to that of an airfoil executing vertical transla- 
tory oscillations, with downward displacement 


bA = ball 
Oscillatory lift and moment are given by 


K® = xpU°b(ak,)H 


M, = xrpU*b?(am,)H ae 


Finally, by combining Eqs. (A13) and (A14), the 
following result is obtained: 


K = K + K® = apU°%bk,H 
My = Mo? + M,© = 2oU2b2m,H 


where 


kn = alka + 2ik/V1 — M2) 
My, = Qing { 


(A15) 


APPENDIX B.—MUTUAL INDEPENDENCE OF RESONANT 
AND ANTIRESONANT MODES 


Under certain assumptions, based upon considera- 
tions of possible system mode shapes outlined in Ap- 
pendix B, it will be shown that propeller flutter may 
occur in either of two distinct forms, analogous to 
symmetrical and antisymmetrical wing flutter. For 
simplicity the development will be confined to a pro- 
peller with three blades, although the methods are 
general and similar conclusions may be obtained for a 
propeller with any number of blades. The following 
assumptions are made: 

(a) The vibratory torque reaction due to hub inertia 
and effective inertia of the drive system may be neg- 
lected. 

(b) The hub its rigid. 

(c) The propeller shaft is prevented from executing 
lateral or axial motions. 

Displacement shapes for the entire propeller are ob- 
tained by combining modes of the individual blades so 
that boundary conditions are satisfied for the propeller 
asa whole. It is necessary to introduce a complete set 
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—— 
of propeller shapes for each blade mode (consisting of gi = generalized coordinate corresponding _ to 
the maximum number of linearly independent shapes the first resonant flexural mode (zero fre- 
that can be obtained without violating the boundary quency) 
conditions). For present purposes we consider pro- ge = generalized coordinate corresponding to the 
peller displacement forms that are composed of the second resonant mode 
first (zero tenet a — wren modes, the g; = generalized coordinate corresponding to the 
first antiresonant f exural mode, and the first torsional first antiresonant mode 
mode of the individual blades. - ; . : 
gs = generalized coordinate corresponding to the 


From assumption (b) it follows that the angular dis- 
placements about the propeller axis must be the same 
for all blades; only one propeller displacement form 
can be obtained from each of the resonant modes (Fig. 
Bl). 


antiresonant mode (Fig. B2), with one blade bending in 


Two independent forms are obtained from the 


opposition to one or both of the others. 

A maximal set of three linearly independent forms 
can be obtained from the torsional mode; it is conven- 
ient to employ the set shown in Fig. B3. 

As the first step in formulating the flutter equations 
for the entire propeller, we consider the generalized 
forces acting on a single blade with four degrees of free- 


dom; the generalized coordinates are: 


first torsional mode 


The displacement forms associated with these general 
ized coordinates are normalized so that the tip dis 
placement (either linear or angular) is equal to unity 
when g; = 1. Then the generalized forces (including 
inertia, elastic, centrifugal, and aerodynamic effects) 


may be expressed as follows: 


4 
0.= 0 Ag, (BI 
Jj 


1 


The coefficients A;; are complex functions of the rota- 
tion speed, advance ratio, air density, and frequency of 


blade vibration. 


Generalized coordinates for the entire propeller, corresponding to the seven basic shapes indicated in Figs. Bl, 


B2, and B3 are denoted by qu, gz, . . ., gz, respectively. 


If the blades are numbered in the counterclockwise sense, 


starting from the two o'clock position, then the column of generalized forces acting on the jth blade is represented 


by a matrix equation 


Oo 


where 


= Aq (B2) 
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[An Ai (Ai3/2) Ais Ais (Ays/2) Aus 
\ 1») Ax» (Ao; 2) Ao; Any (Avg 2) Any 


(Az3:/2) (Az2/2) (Az3/4) (A3s/2) (Aas/2) (Asa/4) (Asa/2) 





























| A” = A3 A32 (Ag3/2) Ass Ass (Ags/2) Assy 
Ay Ag (Ag3/2) Age Au (Ags/2) Aas 
(A 41 2) (A 49 2) (A 43 ‘4 ) (A 43 2) (A 44 2) (A 14 4) (A 14 2) 
L Ay Ag (Aq3/2) Aa: As (Ags/2) Aas wl 
) J Ay Ax (A 43/2) — Aj; 114 (A 14/2) -Aj, 7] 
A», Ao (Ao3/2) —A2; Ax (A 24/2) —A 
(Aa 2) (Az. 2) (A 33 4) (— 1 22 2) (As 2) ( 4 34 }) ( —Axs 2) 
A o = —A3, —Azxz: (—A33/2) A33 —A3, (—A 34/2) 154 
A4 A (A 43/2) —A4; A 4s (A 43/2) 144 
(A 4/2) (A q/2) (A 43/4) (—A4:2/2) (A 44/2) (A 44/4) ( u/2 
tL —Ay —Ay (—A 43/2) A 4; —A4s (—A44/2) A 4; a 
An Ap — Ais 0 Aj; ~Ajs 07 
A» Ao — Ao; 0 Ao, Avy 0 
—Ax —Az A33 0 —Ax3s A 34 0 
j* = 0 0 0 0 0 0 
————— Ay Aw —Ax43 0 Au —A 44 0 
—A 41 —A 42 Aas 0 —A 44 A 4 0 
0 0 0 0 0 Ou 
to 
_ The flutter equations are 
the Aq =0 (B3) 
the where 
the 3A 3A jy 0 0 3A is 0 0 7] 
342 3Axn O 0 3Ang =O 0 
0 0 3A:3 2 0 0 3A 34 2 0 
eral 4 = A"! A = i QO 0 Q 2A33 O 0 2A 34 
dis 344 3Ay O 0 3Ay 0 0 
nity 0 0 3A43/2 O 0 3Aq4/2 0 
ding LO 0 0 2A,, O 0 2A ” 
ects) 
Three independent sets of equations are obtained from Eq. (B3) 
(B1) 
Angi + A wge + Aisgs = 0, Angi + A22qG2 + Aongs = 0, Aagi + Avge Ausgs = 0 (B4) 
‘ota- 
y of A335 + A; 196 = 0, A 1242 + A 4496 = ( ( B5) 
A33q1 + Assgz = 0, Aags + Augz = 0 (B6) 
It will be observed that Eqs. (B5) and (B6) are equivalent. The implications are (see Figs. B1, B2, B3 and note 
Bl, that 1 and 2 are coupled with 5, 3 with 6, and 4 with 7): 
nse, 
ated (1) There is no coupling between resonant and antiresonant modes. 
(2) The flutter equations may be derived for a single blade. 
(B2) 


(3) It is permissible to investigate separately the possibilities of flutter involving resonant and antiresonant 


modes. 
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The Drag Problem at High Supersonic Speeds 


ADOLF BUSEMANN* 


National Advisory Committee for Aeronautics 


SUMMARY 

The problem is to find the real critical Mach Number in the 
theoretical sense—the first appearance of shock waves in a non- 
viscous fluid following given boundary conditions. It is known 
that some solutions of shock-free potential flow with a mixed 
subsonic-supersonic velocity pattern exist. This fact suggests 
that exceeding the sonic velocity is not critical. It can be shown, 
however, that the genuine potential flows with a mixed pattern 
are the exceptions; the flow past a given body at high subsonic 
speeds leads to a unique solution containing shock waves in the 
overwhelming majority of the cases. This result settles many 
previously unanswered questions in the transonic speed range. 
The main conclusions can be extended to three-dimensional cases. 


INTRODUCTION 


7 SUBJECT OF THIS PAPER is the drag of a single 
body moving with constant speed in an unlimited 
nonviscous fluid.'_ Complications caused by the well- 
known induced drag due to lift are avoided by the choice 
of a symmetrical body at zero angle of attack and 
at zero yaw. The fluid is considered compressible 
and—in order to get a proper classification of the case 
in consideration among all possible cases—the speed of 
the body may be unlimited. There is only one special 
value of the body speed where the steady flow relative 
to the body does not exist in the usual meaning of the 
term; this is the case where the body speed is equal to 
the free-stream speed of sound. If this special value is 
omitted or replaced by a limiting process, the value 
being approached from either side, the speed of the 
body has the simple alternative to be either ‘‘subsonic’”’ 
or ‘‘supersonic.’"’ The perfect flow relative to the body 
in the steady state approaches in its velocity distribu- 
tion the value of the body speed at great distances from 
the body, but the positive and negative excess velocities 
in the neighborhood of the body do not respect any 
limitations connected with the speed of sound. In 
the so-called transonic speed range of the body the 
relative flow is of the mixed subsonic and supersonic 
type. By combining the body speed and the flow type 
the proper classification has the following designations: 
pure subsonic flow, mixed subsonic flow, mixed super- 
sonic flow, and pure supersonic flow. 

This classification does not say that one and the 
same body generally has all four speed ranges; on the 
contrary, a sharp or pointed nose, the requisite for a 
pure supersonic flow, is a severe handicap for a pure 
subsonic flow unless symmetry is used as a pure mathe- 
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matical fact rather than as a dependable physical 
quality. It may be worth while to add that this 
difficulty in the classification and the obvious mathe- 
matical reasons for preferring the two pure flows for 
simplicity are only indirectly connected with the prac- 
tical drag problem. It is true that the pure supersonic 
flow requires a sharp nose and that the blunt nose suits 
the nature of the subsonic flow. The question whether 
and where the shape of an optimum body must change 
from blunt to sharp noses belongs to the mixed flow 
problems and often depends upon the secondary condi- 
tions applied to the desired optimum body. 

One more remark about exceptional cases may be 
permitted concerning a special gas whose derivative of 
pressure with respect to volume for isentropic changes 
of state is a (negative) constant. This gas, called 
“‘Hooke’s gas’”’ or perfect gas with the specific heat 
ratio y = —1, has no mixed flow pattern or transonic 
speed range. It is therefore of no interest whatsoever 
for settling any question not related to the pure subsonic 
and the pure supersonic flow. The problems discussed 
in this paper are solved for gases with concave isen- 
tropic lines throughout in pressure volume coordinates. 
A convex line would reverse some of the alternatives 
(expansion shocks instead of compression shocks), and 
Hooke’s gas must not be considered. 


THE SUPERSONIC ENCLOSURES 


It is more convenient for the discussion to reduce the 
number of variables and to investigate the two-dimen- 
sional flow. In the terms of the two-dimensional flow 
the borderline between a subsonic and a supersonic 
region is called a “‘sonic line,” representing the locus 
where the flow velocity equals the local speed of sound 
of the fluid. In a potential flow this coincidence is only 
possible at a certain state of the gas; therefore, pressure, 
density, temperature, and speed of sound are constant 
along the sonic line and are in the so-called critical or 
sonic state. The body is considered as having a general 
shape, which is independently and arbitrarily given. 
The compatibility of the subsonic and supersonic re- 
gions settles the shape of the sonic line. The sonic line, 
therefore, is the last thing indeed that can be considered 
as known in advance or can be used as a substitute for 
the body shape. There are, however, some general 
rules for tracing the sonic line in a potential flow. 
The key to these rules is the exceptional flow behavior at 
a point (Fig. 1) where the sonic line crosses a streamline 
just at a 90° angle. In the neighborhood of such a 
straight crossing the sonic line has a finite radius of 
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curvature, which depends upon the pressure gradient 
normal to it, whereas the streamlines are closely parallel. 
The finite curvature having the subsonic region always 
on the inner side makes the progress to the sonic line 
change definitely from an oblique upstream direction 
to an oblique downstream direction or vice versa. 
The potential lines orthogonal to the streamlines which 
allow the measurement of the upstream or downstream 
progress will have either no intersection or two inter 
sections with the sonic line; for those potential lines 
that have 
between these two intersections lies necessarily in the 
subsonic region. Considerations of this type add up to 
the general rule: Each potential line can only have 


two intersections, the internal segment 


zero, one, or two intersections with one and the same 
sonic line; if there are two intersections, the internal 
segment between these two intersections must be a 
subsonic one. 

There are two types of potential lines in the mixed 
subsonic flow past a single body (Fig. 2). One type 
passing in front or in rear of the body connects the two 
subsonic infinities; the other type connects a point of 
the body contour with either subsonic infinity. This 
situation excludes internal subsonic segments and 
therefore, the maximum number of 


reduces, inter- 


1949 


sections to one. Only those potential lines that reach 
the body can, have one intersection; the others have no 
intersection at all. The general aspect is, therefore, 
that the supersonic enclosures are always adjacent to 
the body and imbedded by the subsonic flow; the sonic 
lines intersect their potential lines only once. The only 
remarkable variation occurs when the expected super- 
sonic region turns out to consist of a couple of successive 
regions, every one adjacent to the body contour. 


THE CURVATURE OF THE BoDY CONTOUR 


At the sonic line, the two different flow types 
supersonic flow—determine 


The subsonic type is 


the subsonic and the 
the compatibility conditions. 
similar to the incompressible flow, and the supersonic 
type is known by the method of characteristics or Mach 
lines. 
Every two-dimensional supersonic flow gets its 
structure by two families of Mach waves crossing the 
streamlines at the so-called Mach angle, which, in 
potential flow, varies with the pressure or the flow 
speed only. The two families, therefore, cross one 
another at twice the Mach angle, and this internal 
corresponding to Mach’s cone—is bisected by 


the external angle, which is supple- 


angle 
the streamlines; 
mentary to the internal angle, is bisected by the poten- 
The supersonic potential flow is 


‘ 


tial lines (Fig. 3). 
represented by means of the Mach lines, which are not 
only used to show the structure but are also utilized 
as a quantitative system. A selected system of lines 
in both families, for instance, indicates the deflection 
of the streamlines in steps of turning angles of 1°; the 
same system of lines indicates the pressure variations 
by steps up or down in a certain function of the real 
pressure. As velocity direction and pressure are not 
mutually related, the coordination of pressure steps 
and turning angles is opposite in the two different 
families. The well-known roughness of a supersonic 
flow contrary to a subsonic flow makes it necessary to 
a series of compression and expansion lines 
The most dangerous 


expect 
alternating in both the families. 
items with respect to the existence of a potential flow 
are the compression lines, which have the quality 
to collide and to form shock waves. In the conven- 
tional system, the pressure steps upward—or the com- 


pression lines—are drawn as the solid lines—danger 
and the expansion lines are drawn as_ broken 
lines. With respect to the velocity angle, however, 
the broken line in one family corresponds to the solid 


line in the other family, according to the interchange of 


lines 


pressure and angle coordination. 

The normal aspect of a supersonic flow is, therefore, 
that both families have both kinds of lines. If one 
family disappears or has zero intensity, the lines of 
constant pressure are the Mach the only 
existing family (Prandtl-Meyer expansion). If the 
two families exist with opposite signs, the lines of 
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constant pressure intersect the internal angle between 
the Mach lines; if the families have the same sign 
both compressions or both expansions—the lines of 
constant pressure intersect the external angle between 
the Mach lines. 

The sonic line is such a line of constant pressure, but 
it is an exceptional one. On one hand, the supersonic 
region is on one side only, and, on the other hand, the 


Mach angle is 6 = 90°; the internal angle between 
the two families is 28 = 180°, and the external angle is 
180°—28 = 0. Thus, the sonic line can only intersect 


the internal angle, and, according to the previous 
results, the two families of Mach waves must be of 
opposite sign. Moreover, the Mach waves exist only 
on one side of the sonic line, and, consequently, there is 
no choice which family gets which sign. The incoming 
Mach waves must always be expansions and the out- 
going waves must be compressions. This is the one and 
only possibility left for all Mach waves reaching the 
sonic line. When the sonic pressure corresponds to a 
full step in the pressure scale, both families must join 
at the sonic pressure as their final point. As this is an 
extremely suitable convention, the general rule reads: 
At the sonic line of a potential flow, expansion lines 
only can arrive and they fictitiously reflect as com- 
pressions. This structure of a supersonic enclosure is 
demonstrated in Fig. 4. On the wall of the body the 
greater number of expansions in the beginning makes 
the pressure drop until the incoming compressions lines 
become equally numerous. The curvature of the body, 
however, is given by a convex deflection in both the 
outgoing and the incoming family. The opposite sign 
of both families in the monotonic pressure changes 
produces a completely monotonic character for the 
streamline curvature. The supersonic regions can 
exclusively border convex parts of the body. 

The single two-dimensional body cannot be concave 
throughout; nevertheless, a general shape of the body 
has not been restricted hitherto to a convex shape 
throughout. The existence of a potential flow past a 
body of arbitrary shape seems already questionable 
when the body shape is restricted to convex contours. 
The possibility of splitting the supersonic region into 
several successive regions that border 
only the convex arcs prevents an immediate decision. 
It pays, however, to look closely into the existence 


supersonic 


problem. 


INFINITESIMAL VARIATIONS OF THE Bopy CONTOUR 


There is a strong temptation to take a given solution 
of a convex body and to buckle its surface gradually 
until the valleys of the waviness are concave. The 
nonlinear differential equation of the flow indicates 
however, that it is advisable to start with infinitesimal 
variations that are simpler and can be superimposed. 

When the walls of a given incompressible flow are al- 
tered by a local depression or protuberance, the disturb- 
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ance pattern surrounds the source of the trouble and its 
excess velocities decrease with the distance from the 
source. This pattern holds approximately in the 
subsonic flow, the main difference is that the lateral 
disturbance pattern gains weight relative to the longi- 
tudinal pattern with increasing Mach Number. At 
supersonic speeds the whole disturbance pattern is 
restricted exclusively to the internal angle of the Mach 
waves downstream of the* disturbance source. The 
distribution within this angle in two-dimensional flow 
is not at all uniform. A disturbance source inside the 
flow uses its two outgoing Mach waves and a source at 
the wall uses the single outgoing Mach wave as the main 
propagation line of the disturbance. When the small 
depression of the walls has a small finite length, the 
alley between the two Mach lines starting at the two 
ends is the track of the main disturbance, and the dis- 
turbance remains of about the same intensity inde- 
pendent of the distance from the source. The second 
family of Mach waves crossing this alley gets a weak 
disturbance like the shift of a light beam passing 
through plane parallel glass. 

For the following discussion of the mixed flow a 
symmetrical body with a symmetrical mixed-flow 
pattern will be chosen, and the upper surface may have 
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Fic. 5. The disturbed mixed flow past a symmetrical body. 


a single symmetrical supersonic region (Fig. 5). Nowa 
small depression of symmetrical shape may be applied 
at the axis of symmetry where the highest supersonic 
Mach number occurs. When this depression is created 
suddenly, the nonsteady flow within the supersonic 
region also has a restricted disturbance area, and the 
main disturbance travels through the outgoing alley of 
Mach waves. This alley ends on the sonic line and gives 
access through the subsonic region toevery other point in 
the flow. It, therefore, is not absurd to expect finally a 
symmetrical flow pattern. The supersonic flow character 
valid at the disturbance source and its neighborhood, 
however, introduces a severe difficulty because the in- 
fluence of any action is strictly downstream; no coopera- 
tion, suction, nothing can possibly facilitate to have the 
corresponding share of the disturbance pattern traveling 
on the incoming alley. When the symmetrical con- 
tribution happens to come down the incoming alley and 
to take over exactly 50 per cent of the disturbance, this 
phenomenon regarded from the disturbance source 
must be called purely accidental, and, moreover, it 
represents an amazing degree of precision. The way 
across the subsonic portion between the end of the 
outgoing alley and the entrance of the incoming alley 
is the only connection for controlling the exact alley 
and the character of the disturbance sent down on it 
and for checking the results. It is obvious that here a 
symmetry is expected of a process that is not sym- 
metrical in itself. If, on the other hand, the idea of 
symmetry is abandoned, the potential flow loses one 
of its most valuable quantities—namely, uniqueness. 
What may be expected when the normal character of 
subsonic and supersonic flow is considered? The 
supersonic region sends the main disturbance down the 
outgoing alley. The subsonic flow spreads around the 
end of the alley a disturbance distribution with de- 
creasing amplitudes, but a more concentrated reaction 
must be expected to travel down the supersonic alley 
that starts at this end of the first alley (Fig. 6). The 
second alley gets contact to a third alley by the re- 
flection on the solid walls of the body, and a zig-zag 
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disturbance in the supersonic region will be the result. 
It may be an interesting problem to determine the condi- 
tions near the sonic line where the disturbance reverses, 
but the only thing that really matters is the alternative, 
whether the disturbance on this zig-zag decreases or 
increases. 

The first guess may be that the intensity will de- 
crease because of the contributions to the subsonic 
region. If this guess holds true, the original depression 
can be made small enough to avoid shock formation 
before the disturbance dies completely out. In this 
case a potential flow is guaranteed. The potential 
flow, however, is reversible, and there must be another 
possible solution where the disturbance is mainly 
upstream of the depression. As several infinitesimal 
variations can be superimposed, the original flow can 
be superimposed at the same time by a depression with 
an upstream disturbance and by an equivalent pro 
tuberance having a downstream disturbance pattern. 
Both variations combined do change the potential 
flow but do not change the body. The potential flow 
past the body is, therefore, open to a great number of 
variations according to different shapes of the super- 
imposed positive and negative depressions. The 
potential flow of the mixed kind is, in this case, existent 
but not unique. 

The opposite guess—that the disturbance is increas- 
ing in such a manner that, finally, a shock wave must 
occur—has a completely different aspect. Even the 
smallest depression can only lead to a flow that is 
regular everywhere except at the point where the sonic 
line reaches the body contour. This is the point where 
the zig-zag repetitions are infinitely crowded. The 
punctuation of all points where the body contour and 
the original sonic lines intersect, would make a genuine 
potential flow possible exactly in the same manner as in 
the case of decreasing amplitudes of the disturbance 
(Fig. 7). The previous result of existence, but not 
uniqueness, of the potential flow can therefore be trans- 
ferred to the punctuated flow. This transfer leads to the 
chance to settle the problem of selecting the desired 
solution from all the solutions possible. It is obvious 
that smooth or analytic flow at all the punctuated 
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Fic. 7. The punctuated potential flow. 


points means a genuine potential flow. For a given 
body shape different from the original, however, the 
freedom of the nonunique solutions permits a smooth 
flow at either the entrance or the exit corner or a sym- 
metrical flow past a symmetrical body, but it does not 
permit smooth flow at both corners at the same time. 
There is, however, another way to find new potential 
flows with smooth flow on all punctuated points: 
Every time the disturbance is reflected on the body 
contour, a chance is given to counteract the disturbance 
by a suitable shape and to make it disappear. This 
shows that a genuine potential flow exists for the origi- 
nal body and those variations that are self-compen- 
sating. The guess that the amplitudes increase leads, 
therefore, to a restricted number of existing genuine 
potential flows; it keeps the uniqueness and makes the 


existence doubtful. 


THE DECELERATED SUPERSONIC CHANNEL FLOW 


There is nothing like thermodynamics to guide the 
engineer concerning what can and cannot be done by 
The simple rule is to keep the entropy 
sources small. In some idealized processes, such as 
the two-dimensional perfect fluid in pure subsonics, 
no entropy augmentation is possible, and the result 


his own skill. 


is no drag regardless of shape, the well-known “Paradox 
of d'Alembert.” But in supersonics the wave drag 
exists, and this drag really depends on the skill of the 
designer to keep it small. The drag is not entropy 
augmentation right away. The flow at first is covered 
with waves, and, if this pattern is sufficiently long, the 
compressions coalesce and produce shock waves. As 
long as the waves exist, there is still a chance to take 
them out, and every time this is done the drag is di- 
minished. The parallelism of the supersonic flow in 
wind tunnels depends on the designer’s skill and is not 
a consequence of parallel channel walls. Fig. 8 shows a 
parallel supersonic flow in a channel, which is dis- 
turbed by a depression on the wall. The flow reaction 
on the wall is like a drag, while the flow loses momentum. 
When a compensating protuberance takes the whole 
disturbance out, the previously existing momentum is 
re-established and the reaction on the wall is like a 
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thrust. (The compensating deformation of the wall 
is shown as contrary to the original disturbing shape. 
If the distance between the original disturbance and 
the compensating deformation is not a whole number of 
wave lengths, the compensating slope down may 
come first and the slope up later; the compensating 
shape may thus be a repetition of the disturbing shape. ) 
The only important knowledge is that the compensating 
deformation requires a certain length. A _ self-com- 
pensating variation must be longer than this length. 
In the two-dimensional supersonic channel flow the 
different modes are harmonic and have commensurable 
wave lengths, but that is not important. If the walls 
are not carefully shaped or if the flow comes out of a 
nozzle it always takes a certain compensating length 
to get the flow smooth, and the shape within this 
length has to be accurately compensating both in two 
dimensional and axially symmetrical flow. 

The precision of the Laval nozzle to produce a 
parallel supersonic flow is relatively small compared 
with the precision that would be required for the 
reversed flow, the decelerated Laval-diffusor flow. 
Fig. 9 shows a decelerated zone between two parallel 
channels. The homogenous deceleration between two 
straight convergent walls is connected with the two 
parallel channels by a transit are of concave curvature 
at the entrance and a convex are at the exit. The 
first arc produces the compression indicated by the 
thin compression lines; the second are cancels these. 
When such a decelerating nozzle is passed by a dis- 
turbed flow, the disturbance amplitudes are increasing. 
This process is demonstrated in Fig. 10. A depression 
ahead of the deceleration produces a disturbance, and a 
compensating deformation of the walls takes it out 
after the deceleration is undergone. Now, the thrust 
gained by the compensation is much greater than the 
drag of the first depression. The single compression 
line of the disturbance makes the surrounding lines 


more convergent. Thus, the compression line absorbs 
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on its way through the deceleration some of the thin 
compressions and, therefore, comes out of the decelera- 
tion much stronger than before. The single expansion 
line of the disturbance correspondingly makes the 
surrounding lines more divergent and, therefore, rare- 
fies the weak compression on both of its sides (the ex- 
pansion itself is, in reality, a divergent wedge and is 
merely simplified in the figure when represented by a 
single line). The rarefaction of the compression lines 
has the effect that the convex transit are prepared for a 
homogenous arrangement of compression lines pro- 
duces an outgoing expansion line for every missing 
incoming compression line. Thus the intensity of the ex- 
pansion grows also on the way through the deceleration. 

Another way to explain the growing of the disturb- 
ance in a deceleration, which leads to the same results, 
is given by the law of momentum in its longitudinal 
component. The walls of the parallel flow in Fig. 8 
contribution to the tangential momentum 
component. It was, therefore, without any significance 
that the average pressure of the 
different from the undisturbed static pressure on the 
On walls with slope, however, an average excess 
Regardless 


had no 


disturbed zone is 


walls. 
pressure contributes to the momentum. 
of the disturbance shape, the excess pressure is positive 
according to the second order term in a power series 
development for Ackeret’s supersonic relation between 
pressure and local angle of attack. This term is known 
to be always positive for reasonable perfect gases, and it 
is positive, too, at Mach Numbers near 1.0 for all gases 
with a concave pressure-volume plot. This positive 
excess pressure makes the momentum losses increase 
in a decelerated supersonic flow because of its conver- 
gent stream lines. The thrust acting on a compensat- 
ing deformation of the walls is, 
is undergone than 


therefore, 
when a deceleration the 
sponding drag on the first deformation (Fig. 10). 


greater 
corre- 


THE MIXED FLOw IN A CURVED CHANNEL 


The curved channel flow between concentric circles 
has a potential solution, all stream lines being con- 
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centric circles. The centrifugal forces allow a mixed or 
transonic flow with the subsonic region on the outer 
wall. The infinitesimal cyclic disturbances leading to 
corrugated steam lines can be discussed by an additional 
potential, which has an amplitude function in the radial 
direction given by Fig. 11. (This is G. I. Taylor's 
problem of 1930.°) Every maximum or minimum of the 
amplitude function represents an uncorrugated stream 
line which can be used as a solid wall. Every zero of 
this function gives a corrugated stream line of constant 
pressure and can be used for a wall with a thick bound- 
ary layer. In either case, the supersonic region has many 
stream lines available to satisfy the given boundary 
condition, but, in the subsonic region, only one single 
stream line at a time can be made either undeformed or 
isobaric. In Fig. 11 the stream line at the subsonic 
infinity is chosen to satisfy exactly both conditions, 
but practically every stream line outside a minimum 
distance from the sonic line is already extremely close 
to these conditions because of the asymptotic character 
of the amplitude function. This result proves the 
suggestion previously made that the position of the sonic 
line is more significant for the determination of the 
disturbance modes than the proper subsonic boundary 
conditions on the wall. The modes are not harmonic, 
and the sonic line is neither a node nor a loop but close 
to the end of interesting amplitudes. 

The decelerated flow in a curved channel without and 
with superimposed steady disturbances, treated by G. 
Guderley, 1945,' is the problem that really discloses 
the answer to the question of increasing or decreasing 
amplitudes. That his fundamental result came as a 
verification rather than as a surprising discovery, 
because the indications were already strong for the in- 
creasing amplitudes, may serve as the welcome reason 
to skip the details of this troublesome investigation, 
since its application to a straight channel (the actual 
proof of its reliability) is already given. The deceler- 
ated mixed flow that has a fading supersonic region 
must still keep the main disturbance inside this super- 
sonic region. The sonic line has a discharge, and the 
stream lines escaping through it cannot be troubled 
They correspondingly have their 
It may be 


by shock waves. 
momentum losses regained on the way out. 
considered as an open question whether the total losses 
increase or decrease by such a concentration of the dis- 
turbance; in any event, the realization of the losses 
by unreversible shock waves comes sooner. 

The accelerated flow decreases accidental disturb- 


ances; the decelerated flow increases the amplitudes. 
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This reason is sufficient to establish the following 
selection rules for the perfect flow: The entrance corner 
of supersonic regions must be smooth; the exit corner 
may be the seat of a shock wave pattern or of a single 
This rule holds trpe for three-dimensional 
flow because the the ability to 
oscillate and the positive excess pressure on divergent 


shock wave. 
physical factors 
stream lines——are just the same. 

An effect well known to everybody who played with a 
string may be used as a physical analogy in order to 
underline the conclusions. In Fig. 12 a string hanging 
from a hole represents a system able to pendulate in 
the gravitation field. The amplitudes of all its un- 
harmonic modes are decreased by lengthening and 
increased by shortening the string. It is a_ well- 
known effect that the string gets extremely small ampli- 
tudes when the lengthening starts with zero length at 
the hole and proceeds, but when the string is pulled 
back through the hole it nearly always gets extremely 
large amplitudes and ends the movement in a clash. 
The energy is due to pulling against the positive centrif- 
ugal forces of the modes. 

The schlieren photograph (Fig. 13) demonstrates the 
practical case at a Mach Number 0.75 on the NACA 
662-015 airfoil in a Langley wind tunnel. The two 
dotted sonic lines are plotted according to the experi- 
mental pressure distribution gained by pressure holes 
in the wind-tunnel walls. The contrast between the 
flow character at the entrance and the exit corners 
is obvious. One possible objection may be that some 
of the disturbances caught by the spark light are not 
steady but traveling disturbances. The unsteady flow 
character, however, can also be used to explain the 
smoothening effect at the entrance corners and the 
roughening effect at the exit corners: The entrance 
corners are unstable, and the exit corners are stable 
equilibrium position for traveling waves normal to the 


surface. 


THE New ASPECT OF THE SUBSONIC MIXED FLow 


The result that the disturbances are amplified in a 
decelerated supersonic or mixed flow settles the ques- 
tion about existence and uniqueness of a potential 
flow past a body of general shape. It is the existence, 
not the uniqueness, that is doubtful because the co- 
existence of smooth entrance and exit corners of all 
supersonic enclosures can only be satisfied in rare 
exceptions. A punctuated potential flow is possible in 
the neighborhood of existing cases but has difficulties 
about the uniqueness which require, as a unifying 
condition, that the entrance corners must be smooth. 

Normally, shock waves cannot be excluded at the 
exit corners of supersonic enclosures in mixed flows. 
Similar to the procedure in supersonic nonviscous flow, 
the potential flow must be generalized to the ‘“‘perfect 
flow,’ which degenerates to the potential flow every 





Fic. 13. Schlieren photograph with dotted sonic lines added. 


time no shock waves occur. The existence and unique- 
ness problems have to be settled for this flow type, too. 
The perfect flow inherits the pure subsonic and super- 
sonic solutions and has a better chance in the mixed flow 
to copy the existence and uniqueness of the experimen- 
tal flow. To suggest, however, both existence and 
uniqueness throughout the whole range would be going 
too far. There are already exceptions known for both 
existence and uniqueness. The flow past a body has 
difficulties in the rear by the vacuum line that occurs at 
the Prandtl-Meyer expansion and by the subsonic rear 
pressure problem of a partly separating perfect flow. 
The uniqueness cannot be expected for arrangements 
like the flow through a tube with two throats or like 
supersonic intakes. Nevertheless, the existence and 
uniqueness for the perfect flow past a single, simply 
connected body of general but reasonable shape seems 
to be evident by experimental and mathematical ex- 
perience. 

The real critical velocity of a body where d’Alem- 
bert's Paradox predicting zero drag generally ends is the 
first appearance of the sonic state. As far as the gen- 
uine potential flow penetrates the mixed flow, there are 
still some rare exceptions of zero drag according to 
d’Alembert’s Paradox. It is obvious that the entropy 
augmentation for the neighbor solutions of the excep- 
tions is close to zero and, correspondingly, the drag, too. 
It therefore is an impractical task to draw this sensitive 
borderline between the zero and the positively zero 
drag. The real problem is to determine the values of 
drag for given body shapes within the usual order of 
workshop accuracy, and this practical problem promises 
reasonable solutions. It is, of course, more sensitive 
than the usual subsonic or supersonic flow, but normally 
it is not too sensitive. Supersensitive arrangements, in 
case they may exist, are of no practical value anyway, 
and it does not pay to solve them. Some problems are 
even less sensitive than they are usually feared to be, 
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especially at Mach Number 1.0 or the exact sonic speed 
of the body. Ina recent paper,* G. Guderley solved the 
choking pattern of the free stream at Mach Number 1.0 
far from the body. This result is the final step in con- 
necting both the subsonic and supersonic approach to 
the pressure distribution at Mach Number 1.0, a 
problem that the author started in 1930 in his often 
misinterpreted paper.‘ This solution may re-establish 
the confidencesin wind tunnels near Mach Number 1.0 
with finite test sections. It is the unequaled precision 
of the wind-tunnel flow instead of the infinite test sec- 
tion which is really important. Another result of this 
paper by G. Guderley is the finite drag depending on the 
thickness ratio at and near Mach Number 1.0 even for 
the optimum shape. A first value is a drag coefficient of 
0.08 at 10 per cent thickness and corresponding drag 
values proportional to the 5/3 power of the thickness 
ratio because of the transonic similarity. This settles 
the question whether it is possible to reach a Mach 
Number approaching 1.0 with any genuine potential 
flow without drag; it is definitely impossible. 

When the perfect flow takes over the duties of the 
potential flow, the boundary-layer problems have to be 
combined with the perfect flow. If an infinite gradient 
of pressure rise is considered as a sufficient reason for a 
boundary-layer separation, the critical Mach Number 
with respect to separation would be the first appearance 
of the sonic state. This belief, however, is not true. A 
finite small pressure step is at least faired in the neigh- 
borhood of the wall by the finite subsonic boundary 
layer between. A boundary layer not already in a 
critical state will always suffer pressure steps and usu- 
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ally even the normal shocks of a Mach Number 1.2 or 
1.3 and, more easily, the reflected oblique shocks. When 
the perfect flow pattern at the exit corner is too unde- 
termined for an exact boundary-layer treatment, the 
safe way will be to choose the worst possible conditions 
first and to try afterwafds to press the conditions into 
narrower, dependable limitations if it matters. 

The examples given may illustrate the changes of 
methods which go hand in hand with a serious change in 
the physical phenomenon. In an early state like the 
present one, when solutions cannot be presented but 
previous favorite ideas are designated to be abandoned, 
it is desirable to show how strong the foundation of a 
new theory is and that it really is a guide to applying 
the right tools on the right questions. A contribution to 
this task is the purpose of the present paper. 
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Station Functions and Air Density 
Variations in Flutter Analysis 


MANFRED RAUSCHER* 
Massachusetts Institute of Technology 


‘SUMMARY 


Station Functions have been devised as a means of correlating 
the motions of the masses at different stations in a vibrating sys- 
tem. The tie-in is analogous to that between the ordinates of a 
curve at a series of stations in the development of Simpson’s 
Rule. The use of the Functions combines the conveniences of 
working with discrete masses and of having continuously defined 
deflections. The Functions are shown to serve well both for the 
calculations of vacuum modes and for the determination of flutter 
modes. 

Air Density Variations have been considered as a means of fix- 
ing the flutter conditions corresponding to a chosen value of the 
flutter parameter k. The attractiveness of the scheme lies in its 
automatic coverage of the flutter performance at various altitudes 
and in its avoidance of calculations for hypothetical k’s at air 
densities fixed in advance. The method becomes especially valu- 
able when allowance is made for aerodynamic span effects, which 


greatly increase the labor of trial variations of k. 


INTRODUCTION 


y pve TWO STUDIES HERE REPORTED are phases of a 
broad investigation of techniques of flutter analy- 
sis made for the Bureau of Aeronautics of the Navy De- 
partment [Contract Noa(s)-8790]. More complete 
accounts of the work, with full details of sample com- 
putations, are contained in the reports submitted to 
the Bureau. Companion reports cover applications of 
the Station Function method to the vibration and flut- 
ter analysis of a twin-engined aircraft and the develop- 
ment of methods of making and evaluating vibration 
tests on wings with and without concentrated masses 
at subcritical air speeds. All these reports may be ob- 
tained through the Bureau of Aeronautics. 

Contributions to the two studies covered by the 
present paper were made by G. Zartarian, S. M. Haley, 
G. Fotieo, H. J. Cunningham, P. T. Hsu, and E. Moll¢- 
Christensen. Early survey work on the density varia- 
tion scheme was done by J. Stevens, now with Chance 
Vought Aircraft; in some of the later work on the 
scheme, I. Spielberg, from Wright Field, took part. 
Appreciation is expressed to the Bureau of Aeronautics 
for its permission of the use of the material for this 
paper. 


Presented at the Session on Structures—Dynamic Loads, 
Seventeenth Annual Meeting, I.A.S., New York, January 24-27, 
1949, 

* Associate Professor of Aeronautical Engineering. 


(I) STATION FUNCTION METHOD OF VIBRATION 
ANALYSIS 


(1.1) General 


The replacement of distributed masses by concen- 
trated masses is an artifice often used to simplify the 
dynamic analysis of continuous bodies. In comparison 
with continuous-function methods of the Rayleigh- 
Ritz-Stodola type, the scheme is both crude and labori- 
ous, but it has the advantage of requiring little judg- 
ment and of remaining workable under conditions too 
complex to be handled in other ways. A possibility 
was seen that continuous-function features might be 
grafted on the concentrated-mass plan to increase the 
accuracy and reduce the laboriousness of the plan. 
The investigation of this possibility led to the de- 
velopment of the Station Function scheme here pre- 


sented. 


(1.2) Natural Undamped Vibrations of Continuous Bodies 
on Concentrated-Mass Basis 


When the distributed mass of a section of a body is 
replaced by an equal mass concentrated at the center 
of gravity of the section, two primary inaccuracies are 
introduced into the analysis through the implied sup- 
positions that: 

(1) The resultant of the inertia reactions of the 
elementary masses always passes through the center of 
gravity of these masses. 

(2) A concentrated load produces the same deflec- 
tion as a distributed load of which it is the resultant. 

Both these suppositions result in serious errors unless 
the sections are kept relatively short and, hence, rela- 
tively numerous. By a practical rule,' the number of 
sections should be at least double the number of natural 
frequencies to be found, and a doubling or trebling of 
that number of sections is necessary if the mode shapes, 
as well as the frequencies, are wanted with tolerable ac- 
curacy. The method thus requires the solution of 
characteristic equations of whose roots one-half or more 
are useless. As the method is of interest chiefly for the 
treatment of complicated problems with many roots, 
this overhead of useless roots is a severe burden, and 
its possible avoidance by a modification of the method 
was seen as a primary incentive for such a modifica- 
tion. 
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(1.3) Natural Vibration Analysis Based on Interpolation 

for Inertia Loads Between Reference Stations 

The search for a suitable modification of the method 
was guided by the recognition that the improved 
scheme would have to avoid the inaccuracies of the 
simplifying suppositions (1) and (2) of the concen- 
trated-mass plan mentioned in the last Section. The 
analysis accordingly was put back on a distributed-mass 
basis. At the same time, a continuous definition of the 
accelerations of the mass elements was sought by inter- 
polation between the accelerations at the reference 
This interpolation then became the cardinal 
Through it, the inertia loading 


stations. 
feature of the process. 
on the beam at a circular frequency w could be obtained 
as a continuous function of the deflections of the refer- 
ence stations; from the inertia loading, the deflections 
could be computed in terms of themselves, a number of 
equations equal to the number of reference stations 
being obtained through the requirement that the com- 
puted deflection at each station should be equal to the 
original deflection. The characteristic equation for w? 
is of a degree equal to the number of reference stations 
used. Experience indicates that all its roots are amply 
accurate for the usual purposes if reasonable care is 
taken in the interpolation between the original station 
deflection. 

Because of the crucial importance of that interpola- 
tion, a systematic study was made of different ways of 
providing it. 

The simplest plan naturally was to base the inertia 
loads on a linear variation of the acceleration (and, 
hence, of the deflection) over each interval between 
stations. This plan was tried out on a uniform canti- 
lever beam with two reference stations-——one at the tip 
and one halfway out (Fig. 1). The results showed a 
substantial improvement over those obtained on a con- 
centrated-mass basis with the same reference stations: 
The first frequency was off only 1 per cent (compared 
with 10 per cent on the other basis), while the error in 
the second frequency was just over 5 per cent (com- 
pared with 26 per cent on the other basis). Of the 
mode shapes, the first was excellent and the second 
tolerably good, while the concentrated-mass analysis 
had yielded a fair definition only of the first mode. 
Nevertheless, it seemed worth while to seek still better 
results through a more refined interpolation procedure. 

An obvious way to refine the process was to use inter- 
pe lation by curves tailored to the actual boundary con- 
ditions of the problem, including the requirement of a 
proper tie-in between the slopes, curvatures, and rates 
of change of curvature of the curves on the two sides of 
each reference station. For a cantilever beam with two 
reference stations, as considered before, the interpola- 
tion function wanted in the inner interval was thus one 
corresponding to zero slope and zero deflection at the 
root and to a deflection 2; at Station 1 (Fig. 2), while 
the outer interval required a function that would give 
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the deflection 2, at Station 1 and a deflection 22, with 
zero curvature and zero rate of change of curvature, 
at Station 2. To the seven conditions thus laid down, 
there remained to be added the requirements for con- 
tinuity in the three lowest derivatives of the functions at 
Station 1, making ten conditions in all. These condi- 
tions could have been satisfied in full by a pair of fourth- 
order parabolas. It was decided, however, that atten- 
tion to the third derivatives, both at Station 1 and at 
Station 2, might properly be omitted in functions 
wanted merely for interpolation purposes. The func- 
tions actually tried were therefore third-order parab- 
olas. For a uniform beam, their use yielded, for the 
two frequencies, values that differed only by 0.2 and 
The shape of the 
first mode was given accurately; that of the second 
mode was off perceptibly, but still not seriously (Fig. 5). 
In view of the deliberate violation of some of the bound- 
ary conditions by the interpolation function used, these 
results were considered remarkably good. 


1.2 per cent from the exact ones. 


(1.4) Station Functions 


In the fitting of an interpolation function f(y) to the 
deflections 2; Zz, at a series of reference stations 
V1... Vn, it is always found that the function assumes 


the form 
f(y) = aifi(y) +... + Snfn(y) 


with fi(y) ... fn(y) functions of y determined by the 
boundary conditions but independent of 2 ... 2,. 
Each of the reference deflections z; ... z, thus makes its 
own contribution to f(y); and, as f(y) must satisfy the 
boundary conditions regardless of the values of 2 

Z,, the component functions fi(y) ... fn(y) must indi- 
vidually satisfy those conditions. 

In order that f(y) may have the value 2; at the station 
y1, irrespective of the values of 22 ... Z,, it is evidently 
necessary that fi(yi) = land fo(yi1) = ... = fal) = 0. 
Similarly, it is necessary that each of the other compo- 
nent functions have the value of unity at the station with 
the same subscript as the function while going to zero 
at all the other stations. Each component function 
may thus be viewed as a kind of extension function for 
the station deflection with which it is associated. For 
brevity’s sake, the name ‘Station Function’’ has been 
adopted for such a function. 

The simplest type of Station Function corresponds to 
linear interpolation between station deflections. It is 
represented by a triangle of unit height, with its vertex 
at the station to which the function belongs and with its 
base corners at the next stations to the right and to the 
left. For a station at the end of a range, the triangle 
becomes one-sided. Fig. 3 shows how this scheme, by 
superposition of triangles representing products of sta- 
tion deflections and Station Functions, allows the build- 
ing up of a polygonal curve through any combination 
of arbitrary station deflections. 
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Fic. 5. 


First mode of uniform beam in free bending vibrations. 
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To an interpolation by cubic parabolas, such as was 
considered at the end of the last Section, there corre- 
spond Station Functions that are composed of cubic 
parabolas also. As a result of the continuity of slope 
prescribed at the reference stations, every Station Func- 
tion now runs across the whole range of station inter- 
vals. Fig. 4 shows the graphs of the two functions for 
the two-station setup of the last Section. 

Because the Station Functions are simply components 
of the complete interpolation function, their use in con- 
junction with each other is equivalent to the use of the 
complete single function. The advantage of recogniz- 
ing them individually lies, therefore, not in an improve- 
ment of the results that can be obtained but in a gain 
in the convenience of obtaining those results. This 
gain will be found especially important in the consider- 
ation of flutter with allowance for aerodynamic span 
effects in Section (1.8). 


(1.5) Station Functions Based on Elastic Properties of 

System 

Unlike the simple beam that has been taken for the 
illustrative calculations in the preceding Sections, the 
typical structure encountered in aeronautical practice 
has sharp and erratic variations in its elastic and inertial 
properties from point to point. These variations 
greatly increase the differences between the vibration 
characteristics of different systems, and they make it 
impossible for a set of standard interpolation functions 
to work out equally well in all cases. The attainment 
of a desired accuracy of results then requires either the 
use of a larger number of reference stations than would 
be needed with a system with smooth properties or the 
use of interpolation functions specifically adapted to the 
peculiarities of each system. The benefits of keeping 
the number of reference stations to a minimum are be- 
lieved to make the second alternative the preferable one, 
especially since those benefits carry over in magnified 
importance into flutter calculations. 

A plan for a methodical adaptation of the Station 
Functions to the special features of a system can be 
worked out by reasoning as follows: 

As the system vibrates, the forms it assumes depend 
both on its elastic and on its inertial properties. What- 
ever the latter are, they give rise to a distributed load 
of inertia reactions on the system. Although this load 
will, in general, be of irregular form, the larger features 
of its distribution can be approximated by a polygonal 
curve through suitable ordinates assigned to the curve 
at the different reference stations. If this polygonal 
loading is broken down into triangular loadings, each 
falling from a maximum intensity at one station to zero 
at the neighboring stations, there can be calculated as 
many different deflection curves as there are reference 
stations. These curves, in different superimpositions 
and normalized to give unit deflection at one station at 
a time while giving zero deflection at all other stations, 
then yield a set of Station Functions adapted to the 
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elastic peculiarities of the system. In the functions 
so found for the torsional deformations of a body with 
abrupt changes of torsional stiffness, for example, the 
adaptation will automatically include the proper dis- 
continuities of slope at the points where those changes 
occur. 

To go still further and adapt the Station Functions to 
the inertial, as well as the elastic, properties of the sys- 
tem, a simple plan is to multiply the ordinates of the 
loading triangles by the running mass at every point. 
The loading curves then reflect all the irregular details 
of the actual inertia loading curve, and the deflection 
curves based on these loading curves accordingly reflect 
both the elastic properties and the inertial ones. To 
take care of large concentrated inertia loads, such as 
result from the presence of engines and tip tanks, the 
most convenient plan is to provide supplementary 
functions based on those loads for the stations at which 
the loads are located. For rigidly mounted masses, a 
permanent combination of the concentrated-load and 
distributed-load functions for each station is possible. 
For masses mounted flexibly, the dependence of the ef- 
fective mass on the vibration frequency requires the 
two types of functions to be kept separate. 

For a test of the accuracy attainable by the use of 
Station Functions of the type just described, such func- 
tions were computed for the uniform cantilever beam 
on the two-station basis of the previous calculations. 
The mode shapes derived from these functions are 
shown, together with the exact shapes and the shapes 
previously found, in Figs. 5 and 6. Thesnew mode 
shapes are seen to be virtually exact. The two fre- 
quencies corresponding to them differ only by 0.01 and 
0.15 per cent from the exact frequencies. 


(1.6) Return to Concentrated Station Loadings 


In the pure form in which it has been considered so 
far, the Station Function method attends simultane- 
ously to both the inaccuracies of the concentrated-mass 
scheme pointed out in Section 1.2. Because of the con- 
venience of figuring deflections by the use of influence 
coefficients for the reference stations only, instead of 
influence functions for the whole range of station inter- 
vals, it is of interest to see if the elimination of inaccu- 
racy (1) might do enough good to make worries about 
inaccuracy (2) unnecessary. 

To answer this question, a comparative study was 
made of the results of calculations based (a) on concen- 
trated masses, (b) on distributed masses and distrib- 
uted inertia loads, and (c) on distributed masses and 
concentrated inertia loads [obtained by concentration 
of the inertia loads in a manner analogous to the con- 
centration of masses in procedure (a)]. A_ typical 
example [uniform cantilever beam, two stations for 
procedures (a) and (c), exact solution for procedure 
(b)] showed procedures (a) and (c) to result in errors 
of 10 and 2 per cent, respectively, in the value of the 
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first frequency and errors of 26 and 7 per cent in the 
value of the second frequency. The general conclusion 
was that a concentration of the inertia loads, while 
markedly less harmful than a concentration of the 
masses, results in errors too large to be acceptable. This 
conclusion presumes the concentrations of the inertia 
loads to be made over the full intervals between basic 
reference stations. As influence coefficients will gener- 
ally be available for more stations than are needed for 
the Station Function analysis, there is a possibility of 
working with inertia load concentrations over intervals 
smaller than those between the basic reference stations 
and thus of maintaining a satisfactory accuracy in spite 
of the concentrations. 

In this connection, it is worth recalling that the graphs 
of influence functions are similar in nature to deflection 
curves in that they never have discontinuities in their 
ordinates but may be discontinuous in their slopes, cur- 
vatures, etc. A relatively accurate definition of the 
influence curves is therefore provided by a relatively 
small number of points on the curves, and a knowledge 
of the influence coefficients for a reasonable number of 
stations consequently means, in effect, a knowledge of 
the influence functions over the range containing the 
stations. Instead of concentrating the loads to make 
them correspond to the local influence coefficients, one 
can thus just as simply ‘‘spread’’ the influence coef- 
ficients, by interpolation, to make them correspond to 
spread-out loads. The problem of integrating the 
product of a continuous influence function and a dis- 
tributed-loading function into a station deflection is 
most conveniently handled by a product integraph. 
A summation procedure, utilizing ordinary computing 
machines, is also easy and should be found accurate 
enough with summation intervals of perhaps one- 
quarter the distances between the basic reference sta- 


tions 1... Vn. 


(1.7) Generalization of Station Function Concept 


While the illustrations of the Station 
method in the preceding Sections have stressed the use 
of the functions in connection with bending vibrations, 
it has already been noted in passing that the functions 
have an analogous use in the study of torsional vibra- 


Function 


tions. 

For the treatment of vibrations involving both bend- 
ing and torsion, it is possible to superimpose the Sta- 
tion Functions for the two types of motion. An alter- 
native possibility is to abandon the often difficult or 
unnatural separation of bending and torsion and to deal 
with the combined motion by means of Station Func- 
tions for linear deflections at points on two lines, such 
as the front and rear spars, or leading and trailing edges, 
of a wing. This way of defining the combined motion 
is preferred because it leads directly to the more general 
plan of distributing reference stations at the corners of a 
network of lines on an essentially flat surface and de- 
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Fic. 8. Pyramidal loading for generalized station functions 


scribing the most general weaving motions normal to the 
surface by the linear deflections at the reference sta- 
tions. 

Under this generalized concept, the Station Functions 
become means of defining surfaces, instead of lines as 
before. A suggested way of determining the functions 
is to consider the weaving surface under a succession of 
pyramidal loadings, falling from a maximum value at a 
station to zero on the neighboring mesh lines on the 
four sides of the station (Fig. 8). The number of these 
loadings will be equal to the number of stations. By 
different combinations of the loadings, it is therefore 
possible to find for each station a deflection function 
that gives unit deflection at that Station and zero de- 
flection at every other station. This function is the 
Station Function for the station, based on the elastic 
properties of the system. A function taking into ac- 
count inertia, as well as elasticity, could be found from 
loadings in the form of pyramids multiplied in their 
ordinates by the local mass per unit base area every- 
where. 

An immediate departure from the idea of consider- 
ing deflections in terms of bending and twisting appears 
indicated in the analysis of wings with large structural 
cutouts and sweep. As long as the ribs of such wings 
can be considered rigid in the rib planes, the possible 
deformations of the wing are limited to those of a de- 
velopable surface. The complexities suggested by the 
term ‘‘weaving motions”’ are therefore not yet present. 
In a mild form, they are, however, not far off in large 
flying wings, and they would promise to be prominent 
in ornithopters. Meanwhile, well-developed weaving 
motions already afflict structural components of air- 
craft, such as floors and walls. 


(1.8) Station Functions in Flutter Analysis 


The Station Function method can be used just as 
well for the analysis of flutter as for the analysis of still- 
air vibrations. The same Station Functions serve in 
both cases. Through these functions, flutter calcula- 
tions can be based directly on fundamental elastic, 
inertial, and aerodynamic data without necessary refer- 
ence to parallel calculations that may be made for the 
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still-air modes. The only difference between the prob- 
lems in still air and in moving air is that the absence of 
a wind allows a straightforward solution for the fre- 
quencies and for the mode shapes in the manner seen, 
while the presence of a wind necessitates a cumbersome 
trial-and-error search for a flutter condition, as with 
other methods of analysis. 

It is an inherent weakness of the conventional 
method of defining the deformations of a system by 
isolated station deflections that the method cannot al- 
low for aerodynamic span effects, because these effects 
can only be determined through some kind of interpola- 
tion between the motions at the different stations. 
Such an interpolation is provided by the Station Func- 
tions. Each of these functions gives rise to its own 
spanwise load distribution, and, as the function has a 
fixed form, the load distribution associated with it is 
also fixed. The adjustment of the loads to the resultant 
motion consists simply in the superimposition, in proper 
magnitude and phase relations, of the component load- 
ings determined by the individual Station Func- 
tions. 

The use of Station Functions in this connection works 
out very much like the use of natural mode shapes as 
components of the resultant deformation of the sys- 
tem. It is believed, however, that Station Functions 
will be found to serve better than natural -modes for 
two reasons: 

(1) The introduction of the aerodynamic effects 
does not change the setup of the problem basically from 
what it is in the absence of these effects. A routine 
already established for still-air calculations thus car- 
ries over into the flutter calculations. Work based on 
natural modes does not tie in correspondingly with work 
done in the determination of these modes. 

(2) All Station Functions except those for the tip de- 
flections go to zero at a wing tip. Large aerodynamic 
span corrections at the tip (where these corrections 
give most trouble) are thus confined to the loadings as- 
sociated with the functions for the tip deflections. 
Natural modes all involve tip deflections different from 
zero and so make as much as possible of the aerodynamic 
tip difficulties. If such modes are to be used, it is sug- 
gested that they be put into combinations correspond- 
ing to Station Functions, in the manner illustrated 


by Figs. 4 and 7. 


Because of their ability to describe the deformations 
of a system with great generality, the Station Functions 
appear suitable for the analysis of transient as well as 
periodic motions. They may thus prove helpful in 
work on such dynamic and aeroelastic problems as the 
response of a structure to landing and gust loads and 
the subsidence of displacements produced by measured 
jerks in flutter testing at subcritical speeds. A particu- 
larly simple and effective use of the functions is an- 
ticipated in the treatment of the steady-state problems 
of wing divergence and aileron reversal. 
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(Il) Ark DENSITY AS A VARIABLE IN FLUTTER 
CALCULATIONS 


(2.1) General 


The idea of treating the air density as a disposable 
parameter in flutter calculations was first considered 
as part of a plan of working from vacuum modes to 
flutter modes through stepwise air-density adjustments. 
This plan contemplated the use of Taylor expansions 
of the flutter characteristics in terms of the air density 
at the successive stages and extrapolation from each 
stage to the next. Difficulties with infinite derivatives 
were found to hinder the crucial start from zero density; 
but, once a finite density had been reached, the process, 
while rather cumbersome, proved readily workable.* 
Turning away from the vacuum start idea for the time 
being, the study then went on to explore ways in which 
the density variation plan might be made an economical 
part of regular flutter computations. 


(2.2) Nature of Influence of Air Density on Flutter 


General guidance for the study was sought from an 
analysis of the effects of air-density changes on the 
flutter performance of a typical wing in bending and 
Calculations were made for a model wing with 
Starting from 


torsion. 
the characteristics listed in Table 1. 


TABLE 1 
Wing Characteristics Underlying Illustrative Calculations (in 
Theodorsen’s Notation) 





b = Sin. M = 0.00862 slugs per ft. 

a = —0.30 x = 0.151 (for p = 0.00238) 
Xe = 0.190 w, = 22.17 sec.— 

Toa = 0.636 Wa = 63.43 sec.— 
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sea-level conditions (density ratio p/ Preference = 9 = 1), 


the analysis covered the range of densities down to a 
vacuum. The results of the calculations are shown in 
Fig. 9. 

Worthy of special note are the precipitous drops of the 
curves of the flutter frequency w, and of the flutter 
parameter k at small densities. These drops, together 
with the infinite derivative of the Theodorsen function 
G at the value of k = 0 toward which the k-curve falls, 
account for the difficulties of making the transition to 
flutter from the vacuum vibration. The same drops, 
however, are indirectly helpful for the following rea- 
sons: 

(1) In the case of wy, the rapid final drop to the 
natural coupled frequency w, near 6 = 0 concentrates 
in a relatively small range of low densities most of the 
change in the flutter frequency between sea-level condi- 
tions and a vacuum, so that only gentle variations in 
w, can occur over the range of higher densities from 
standard down to around a fifth of standard, in which 
flutter calculations are commonly made. This insen- 
sitivity of w, to variations in 6 evidently facilitates the 
tracking down of the flutter condition at a new density 
from a known flutter condition at a previous density. 

(2) In the case of k, the drop toward zero at a 
vanishing air density means that a flutter condition 
exists at some density for any & that may be picked be- 
tween zero and a value roughly estimated to be a possi- 
bility at sea level. By taking a & deliberately low at 
the outset and then choosing higher k’s according to the 
results of that first trial, it is thus possible to trace flut- 
ter from some unspecified low initial density up to any 
density that may be of actual interest. 

Both these observations played key parts in the 
working out of the air-density variation schemes de- 


scribed hereafter. 


(2.3) Solution of Flutter Determinant with Air Density 
as Variable 


For flutter of the bending-torsion type, which will 
continue to be taken for illustration, the equations of 
motion can be written 

Az + Bea = 0 

Cz + Dea = 0 
with z and a the displacements in bending and torsion, 
¢, the phase difference between these 
displacements; and A, B, C, D, coeflicients that may 


be expressed in the forms 


A — A, oo A a AQ 


respectively ; 


B = B, + Bo 
c _ C; + C0 


where 0 = p/ Preference, aS before; Q = (w,/wy)?, with we 
the undamped natural frequency of the system in tor- 
. D; are coefficients determined by the 


sion; and A,.. 
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properties of the system, Dy preference and by the flutter 
parameter k. The problem is to find a combination of 
values of 6, 2, and & that will satisfy the requirement 


= AD — BC =0 


The usual procedure is to consider @ as fixed and to 
hunt for a pair of values of 2 and k that fit the pair of 
simultaneous equations given by the real and imaginary 
parts of the expanded determinant. Although the pro- 
cedure is inherently awkward on account of the compli- 
cated manner in which trial changes of & alter the equa- 
tions, the method is the most economical one when the 
flutter condition is wanted only at one particular air 
density. Usually, however, the density of greatest 
flutter hazard is not known in advance, and information 
regarding the variation of the hazard with the density 
would be welcome if it could be had at not too great an 
extra expense. The method of picking a value of k and 
solving the equations for 6 and Q yields that extra in- 
formation for three densities at a cost of about 50 per 
cent of extra labor, while the method of fixing three 
values of # in advance would require extra work amount- 
ing to about 200 per cent. The advantage of the den- 
sity variation scheme is found to lie not in a reduction of 
the number of calculations but in the economy of get- 
ting by with a smaller number of trial values of k. This 
economy greatly increases in relative, as well as 
absolute, importance when three-dimensional instead 
of two-dimensional aerodynamic coefficients are 
used. 

For a simple illustration of the convenience of working 
with a fixed k, bending-torsion flutter with zero struc- 
tural damping will be considered. The absence of 
such damping makes the elastic coefficients A; and D; 
in the flutter expressions real. As the inertia coefti- 
D, are also real, only the aerodynamic 
The expansion 


cients A... 
coefficients A,...D2 remain complex. 
of the determinant thus yields a real equation of the 


form 


R(0, 2) = R; + RO + RQ + Ryo? + RQ + RQ? 


and an imaginary equation of the form 


1(0, 2) = 10 + [02 + 1,02 = Of, + 1,0 + 1,2) = 0 


As the possibility @ = 0 is not admitted in the present . 
discussion, which specifically steers clear of the vacuum 
condition, the imaginary equation gives 


It therefore permits the real equation to be turned into 
a quadratic in @ alone or in 2 alone. The choice of a 
value of k, which fixes the coefficients R,;...Rs and 
I,...Is, then yields immediately the values of @ and 2 
corresponding to that value of k. Only one of the two 
sets of roots for 6 and 2 will usually be found to corre- 
spond to real physical conditions. 
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(2.4) Correlation of Flutter Performance at Different Air 

Densities 

Once a flutter condition has been established at one 
air density not too close to zero, it is generally possible 
to determine the flutter conditions at other densities 
by a convenient process utilizing the information al- 
ready in hand. The basis of this process is the obser- 
vation that the flutter frequency is almost the same 
over a wide range of air densities. 

If the flutter frequency were strictly independent of 
the air density, the indication would be that the aero- 
dynamic forces remained in a fixed relation to the elastic 
and inertia forces at all densities. The aerodynamic 
forces would thus have to be unaffected by changes in 
the air density. The fact is that these forces cannot 
be immune to such changes. For, if the frequency is 
fixed and the air speed increases to make up for a de- 


crease in p, and vice versa, the parameter k = wb/V 


must vary and, hence, alter the relation between the 
real and imaginary components of the Theodorsen 
function C = F + iG, with the result that the air forces 
will change in phase, whatever they may do with regard 
to magnitude. Changes of magnitude are actually in- 
evitable, too, because of the need for the components 
of the flutter motion to adjust their amplitudes, as well 
as their phases, if the work done by the forces over a 
cycle of the motion is to remain zero. The essential 
constancy of the air forces indicated by the constancy 
of w at flutter at different air densities can therefore be 
pictured only as the result of a statistical balancing of 
increases and decreases in the principal air forces. As 
the typical force varies as the product of the dynamic 
pressure pV?/2 and the magnitude (C) = VF? + G 
of the Theodorsen function, the conclusion then is that 
this product at flutter does not vary perceptibly with 
p. With w constant, the quantity (pV?/2)(C)/w*b? = 
(C)p/2k? thus should also come out the same at dif- 
ferent densities. This means that a plot of k?/(C) 
vs. p or @ should be a straight line through the origin of 
the drawing. 

In actuality, w,is not strictly constant over any range 
of p’s, and k?/(C) cannot be expected to be exactly 
proportional to p. But if a density p: differs only mod- 
erately (perhaps by 50 per cent or less). from a density 
p, at which the flutter performance is known, the flutter 
performance at p2 can be found through a linearized 
procedure based on the use of the straight-line approxi- 
mation for k*/(C) for close starting estimates as follows: 

(a) Specify the density pe, or density ratio 6, with 
the understanding that adherence to it will not be rigid. 

(b) From the straight line drawn from the origin 
through the point {6;, [k?/(C)]i} in the @ — [k?/(C)]- 
diagram, read off the approximate value of k?/(C) 
corresponding to #2. This value of k?/(C) fixes the value 
of k for which the new analysis will be carried through. 

(c) At the & thus fixed, express the density ratio as 
# = 6,(1 + 6) and the square of the frequency ratio as 
2 = (1 + €). The flutter determinant can then be 
set up with 6 and ¢ as the only unknowns. 
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(d) In the expectation that 6 and ¢ will prove to be 
small, expand the determinant neglecting squares and 
higher powers of these quantities. 

(e) Solve the pair of real and imaginary equations 
thus obtained for the two unknowns, and, through the 
relations of Step (c), calculate @ and Q. 

(f) If 6 or e, or both, come out too large to let the 
linearization of the process seem justified in retrospect, 
use the values of 2 and @ found according to Step (e), 
and repeat the process with these values in place of 
6 and Q). 

(g) By entering the new point into the {@ — 
[k?/(C)]}-diagram, one can get an idea of the departure 
of that diagram from the assumed straight line for the 
system under study. The new value of © similarly 
indicates how the curve of 2 vs. @ runs between the 
previously known points of Q = 9), at @ = 6, and Q = 
Q,, at @ = 0 (the steep ascent of the curve to the latter 
point is a feature to be remembered in this connection).* 
The two curves are aids toward still closer forecasts of 
the flutter conditions at any new densities for which 
calculations may be wanted. 

For an indication of the closeness with which the 
[k2/(C)]-curve may be expected to approximate a 
straight line in practice, the [k?/(C)]-curve corre- 
sponding to Fig. 9 is shown in Fig. 10. 
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* This feature is peculiar to the case without structural damp- 
ing. When structural damping is present, wy does not tend to- 
ward w, as @ goes to zero. 








Compression Tests of Curved Panels with 
Circular Hole Reinforced with Circular 


Doubler Plates’ 
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ABSTRACT 


A set of 14 panel specimens was tested in compression to de- 
termine the effect of sheet curvature on the strain distribution 
around a circular hole in the sheet. Specimens with no rein- 
forcement around the hole and specimens with doubler plate 
reinforcement of the hole were included. The radius of curva- 
ture of the sheet ranged from 10 in. to infinity. Identical speci- 
mens were tested with and without reinforcement in order to 
evaluate the effectiveness of the reinforcement. 

For the unreinforced specimens, it was found that the plane 
stress theory for flat sheet predicted the mid-thickness strain 
distribution within the errors of observations. 

For the specimens with reinforcements, the plane stress theory 
gave strains that, in the vicinity of the hole, were lower than 
those observed at the extreme fiber of the sheet and higher than 
those observed at the extreme fiber of the reinforcement. This 
was attributed to the fact that the reinforcement was not integral 
with the sheet and did not carry its share of the load. This is 
borne out further by the observation that the shortening of the 
hole diameter was essentially the same for panels with and with- 


out reinforcements. 


INTRODUCTION 


—— PROBLEM of the stress distribution near a rein- 
forced circular hole in a flat plate under stress has 
been solved approximately by Timoshenko,’ who 
treated the reinforcement as a curved beam of constant 
section. Other solutions of this problem based on the 
plane stress theory were obtained by Sezawa and Kubo,’ 
by Gurney,’ and by Beskin.*‘ 

The tests described in this report were carried out at 
the National Bureau of Standards for the Department 
of the Navy, Bureau of Aeronautics, as part of an in- 
vestigation of the effect of reinforcement on the stress 
distribution near holes in monocoque structures. 

Other parts of this investigation are described in 
references 5,6, and 7. Reference 5 gives a plane stress 
solution for the stresses and displacements in the 
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neighborhood of a reinforced circular hole. It gives, 
in addition, the results of tensile tests of a sheet with a 
reinforced hole for which the reinforcement was at- 
tached by rivets. One of the principal conclusions was 
that at least two closely spaced rows of rivets were 
necessary to make the sheet act integrally with the rein- 
forcement. 

The effect of reinforcement in stabilizing the edge 
of the hole and preventing buckling is considered from 
a theoretical point of view in reference 6. It was found 
that, even with large holes, square panels subjected 
to compressive loads are substantially stabilized by a 
reinforcement. 

The stress distribution in the neighborhood of a pin- 
reinforced circular hole was considered in 
reference 7. A plane stress solution was derived and 
found to give stresses that agreed closely with meas- 
ured values in a “‘metalite’’ panel having a reinforce- 
ment bonded to the sheet. 

In aircraft, the sheet in which a circular hole is cut 
is usually curved to a radius that may vary over a wide 
range. The present series of tests was made to in- 
vestigate the effect of curvature of the sheet on the 
strain distribution around a circular hole in a sheet 
under uniform axial load when the hole is unreinforced 
and when it is reinforced with either a symmetrical or 
an asymmetrical doubler plate. The plane stress 
theory, outlined in reference 5, was used to compute the 
stresses and displacements in the neighborhood of the 
hole. The test results are compared with the computed 


loaded 


values. 
SYMBOLS 
A = cross-sectional area far from hole 
P = load 
E = Young’s modulus 
€ = median fiber strain 
€. = median fiber strain far from hole = P/AE 
7. = compressive stress far from hole = P/A 
2uq = shortening of diameter of hole 
r = radial distance from center of hole 
@ = angle between radius vector to (7, 0) and direction of 
load 
v = Poisson’s ratio, 0.3 
a = radius of hole 
6 = outer radius of reinforcement 
t = thickness of panel in unreinforced region 
i i thickness of panel in reinforced region 
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A, Bt, 2.8. 
specimen, determined from Fig. 3 ot reference 5 


constants, depending on dimensions of 


Dimensions of Specimens and Maximum Loads 















Subscripts 
1 = unreinforced region of panel 
j 2 = reinforced region of panel 
| - 
(1) DESCRIPTION OF TESTS 
(1.1) Specimens 
The 14 curved panels tested in this investigation were 
} fabricated at the Naval Air Material Center in Phila- 
delphia. The dimensions of the panels are given in 
2S ' , : ~* T ° ie a 
- [able 1 and in Fig. 1. The panels were of 2458-1 
os aluminum alloy, */;. in. thick, 14 in. long, and a de- 
- veloped width of 24 in. The diameter of the hole was 
” nominally 3in. Six of the panels with no reinforcement 
“ and six of the panels with a doubler plate reinforcement 
= on the concave side of the panel only had radii of curva- 
ture of 10, 20, 40, 80, 120 in., and infinity, and the two 
5 remaining panels with reinforcements on both the con- 
” cave and convex sides had radii of curvature of 10 in. 
id and infinity. The doubler plate reinforcements had a 
‘d volume equal to twice that of the material removed 
® from the hole. Each end of the specimen was held to 
radius by two 1!/4- by 1'/4- by '/4-in. 24S-T aluminum- 
a- alloy angles. 
in : Ténsile tests and single thickness compressive tests 
id were made on specimens from the flat sheet coupons of 
S- the material used in the panels and doubler plate rein- 
E- forcements. ey resulting stress-strain curves are pre- 
sented in Fig. 2, and material properties are given in 
it Table 2. Compressive tests of cylindrically formed 
le coupons showed no significant changes in the properties 
1- of the material. 
le 
st (1.2) Testing of Panels 
d The angle-reinforced ends of the panels were ground 
a flat and at right angles to the axis. The panels were 
“ centered on ground steel blocks (A, Fig. 3) and mounted 
e in a horizontal hydraulic testing machine of 2,300 kips 
e capacity. Plaster of Paris (B, Fig. 3) about '/, in. 
d thick was cast under a small initial load between the 
TABLE 
Radius Sheet Cross- 
of Cur- Thick Sectional Diameter 
vature ness Width Area of Hole 
' Panel (In.) (In.) (In.) (In.?) (In. 
l 10 0.187 24.00 1.474 2.99 
2 20 0). 188 24.00 4.519 3.02 
3 410) 0.187 24.00 4.478 2.96 
t 80 0.186 23.87 4.434 2.96 
5 120 0.186 23.90 $.437 2.96 
f 6 x 0.187 23.98 t. 487 2.97 
10 0.186 24.00 t. 462 3.02 
0.186 23.98 1.468 3.00 
0.187 23.99 4.490 3.05 
0.185 24.03 4.449 3.04 
0.186 24.03 +.479 3.03 
0.185 23.97 4.439 3.04 
0.187 23.86 4.455 3.01 
0.187 24.00 4.496 3.01 
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Nominal dimensions of panels. 
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Stress-strain curves of 24S-T aluminum alloy used in 
Materials A, B, and C are used as 


Fic. 2. 
panels and reinfotcements. 


indicated in Table 2. 


and the steel 
block to take up any lack of parallelism between the 
specimen and the testing machine heads. 

The edges of the panel parallel to the load were 
supported by grooved steel bars C to maintain their 


moveable head of the testing machine 


straightness under load. 


l 


Reinforcement 


Thick- Maximum Maximum 

ness Width Load Stress 

(In.) (In.) (Lb.) (Lb. /In.?) 
None 180,200 40,300 
None 146,400 32,400 
None 117,700 26,300 
None 93,800 21,200 
None 89,700 20,200 
None : cei vibiew ‘ah 
Sym. 0.376 0.653 195,700 43,900 
Sym. 0.369 0.632 79,200 17,700 
Asym. 0.377 0.604 194,000 43,200 
Asym 0.384 0.591 175,600 39,500 
Asym. 0.382 0.611 111,900 25,000 
Asym. 0.381 0.605 103,000 23,200 
Asym. 0.381 0.628 92,300 20,700 
Asym. 0.382 .638 77,600 17,300 
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TABLE 2 
Mechanical Properties of Sheet and Reinforcement : 
—Young’s Modulus— —Yield Strength—-— 
(Kips/In.?) (Kips/In.?) Tensile 
Used in Direction Compres- Compres- Strength Elongation ‘ 
Curve Panels of Rolling sion Tension sion Tension (Kips/In.?) (Per Cent) » 
A 1. 2 oe 2 Long. 10,700 10,500 45.9 aie 73.7 19.0 
8, 9, 10, 14 Trans. 10,700 10,700 50.6 46.4 71.3 19.0 > 
Reinf. on 
panels 7,8 
B 4, 5, Ti, 12, Long. 10,600 10,500 44.5 a 72.2 21.5 
13 Trans. 10,700 10,600 49.5 45.8 69.9 21.5 
% Reinf. on Long. 10,600 12.8 
panels 9, Trans. 10,700 47.7 
10, 11, 12, 
13, 14 
TABLE 3 
Strain Ratios for Panels with Unreinforced Holes 
Gage 
General Length Gage Coordinates* -Radius of Curvature (In.) ——-— ~ 
Location (In.) ¥ y 10 20 40) 80 120 x 
Strain Ratio, e,t/e. 
Unreinf. portion 0.25 1.50 0.00 3.16 3.28 2.80 2.07 2.94 3.15 
0.25 1.75 0.00 2.05 2.03 1.92 Liat 2.03 2.13 
0.25 2.50 0.00 1.20 1.22 Las 1.15 1.29 1.34 
0.50 5.50 0.00 1.03 1.06 1.08 1.12 1.12 1.14 
0.50 9.00 0.00 1.07 1.09 1.09 1.07 0.98 
0.25 0.00 3.50 0.45 0.54 0.42 0.47 0.54 0.50 
0.25 0.00 5.00 0.64 0.67 0.66 0.70 0.70 
0.50 5.50 5.00 1.08 1.07 1.05 1.038 1.08 
0.50 11.00 5.00 ap 0.99 
Strain Ratio, e:t/e« 
Unreinf. portion 0.25 0.00 1.59 —0.85 — 0.87 —0.90 —0.82 —0.44 —0.38 


"2am 


centerline parallel to direction of loading, in. 
€: = strain parallel to transverse centerline. ¢, = 


distance from center of hole along centerline transverse to direction of loading, in. y 


= distance from center of hole along 


= strain parallel to longitudinal centerline 


TABLE 4 
Strain Ratios for Panels with Holes Reinforced by Symmetrical Doubler Plate 
r -- Radius of Curvature (In.)—— a + 
Gage 10 ~ . —— © ——— ~ 
General Length Gage Coordinates* Concave Convex Concave Convex 
Location (in.) x y Sheet reinf. reinf Sheet reinf. reinf. 
Strain Ratio, eyt/e« 
Unreinf. portion 0.25 1.50 0.00 2.38 2.26 
0.50 5.50 0.00 1.06 1.09 
0.50 9.00 0.00 1.02 1.01 
0.25 0.00 3.50 0.54 0.61 
0.50 0.00 5.00 0.79 0.80 
0.50 5.50 5.00 1.02 1.03 
Reinf. portion 0.25 1.60 0.00 1.39 0.71 £23 0.99 
0.25 2.04 0.00 0.00 0.42 —0.07 —0.16 
Strain Ratio, e-t/eo 
Reinf. portion 0.25 0.00 1.60 —0.81 —0.24 —0.43 —0.29 
0.25 0.00 2.04 0.39 0.66 0.36 0.34 








. 
x= 
centerline parallel to direction of loading, in. 
t es = strain parallel to transverse centerline. ey, 


SR-4 wire strain gages '/. and '/, in. long were 
placed on the panel as shown in Fig. 3 for the unreinforced 
panels. The coordinates of the centers of the gages 
relative to the center of the panel are given in Tables 3, 
Two gages were placed on the edge of the 
The gages in the other positions 


4, and 5. 
sheet at the hole. 


were placed back to back, except that, for the panels 
with a single reinforcement, no gage was placed on the 
sheet side opposite the gage near the outer edge of the 


distance from center of hole along centerline transverse to direction of loading, in. y 


= distance from center of hole along 


= strain parallel to longitudinal centerline. 


reinforcement and transverse to the direction of load- 
ing. 

Additional gages were placed 1 in. from each end on 
panel 6 to check uniformity of loading. In the tests 
of all the other panels, the four pairs of strain gages 
5 in. from the transverse centerline and 5.5 in. from the 
longitudinal centerline were used for this purpose. 

Load-strain data obtained were averaged for those 
gages that were placed in symmetrical positions on the 
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Normal deflections of the sheet were measured by dial 
gages D, Fig. 3, mounted on a steel bar E which rode on 
the frame F. The frame was held square against the 
right plate A by leaf springs bearing against the left 
plate A. The frame F was supported by the end rein- 
forcing angles G at the desired point of measurement. 


(2) RESULTS OF TESTS 


(2.1) Strain Ratio 
The median fiber strains were expressed in terms of 
a strain ratio defined as the ratio of the measured 
median fiber strain e, taken as the average strain meas- 
ured on gages back to back, to the average elastic strain, 
= P/AE far from the hole, or 


strain ratio = ¢/e. = ¢/(P/AE) 


The strain-load ratio e/P in the elastic range was 
taken as the slope of the straight line drawn through the 
measured median fiber strain for all gages in a particular 
symmetrically positioned group plotted against load. 
Some bending strain was also present, even at low 
loads. This will be considered later. 

The results of the tests in the elastic range are given 
in Table 3 for the panels with the unreinforced holes, 
in Table 4 for the panels with symmetrical reinforce- 
ment, and in Table 5 for the panels with reinforcement 
on the concave side of the panel only. The maximum 
strain ratio is given in Table 6. It varied from 3.2S 
for panel 2 with no reinforcement and 20-in. radius of 
curvature to 2.26 for panel S with symmetrical rein 


forcement and infinite radius of curvature. 


(2.2) Shortening of Hole and Reinforcements 


The measured shortening of the hole and reinforce- 
ments in the direction of loading was expressed in terms 
of a shortening ratio, defined as the ratio of the meas- 
ured shortening of the diameter of the hole or rein- 
forcement to the computed elastic shortening in the 


absence of a hole, or 


TABLE 6 
Maximum Strain Ratios and Shortening Ratios 


Shortening ratio, 


uaAE/Pa 
Reinf. Maximum Reinforce 
Radius of Thick- Strain ment 
Curvature ness Ratio, Con- Con- 
Panel (in. ) in eAE/P Sheet cave vex 
l 10 3.16 3.11 
°S 20 3.28 
3 40 2 80 2.84 
4 SO 2.44 2.79 
5 120 2.94 3.55 
6 x > 3 15 3 21 . a 
7 10 0.376 2.38 2.92 2.94 2.67 
8 « 0.369 2.26 23.40 2.30 2.66 
9 10 0.377 2.86 3.08 2.42 
10 20 0.384 2.46 2.89 1.61 
11 40 0.382 2.61 2.60 1.84 
12 80 0.381 2.46 3.00 2.07 
13 120 0.381 2.53 3.18 2.60 
14 « 0.382 2.65 3.09 2.00 
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Shortening of diameter of hole, in. 
(b) Panels with 10in. radius of curvature 


Figure 4.- Relation between shortening of diameter of hole 
and load for flat panels and panels with 10 inch 
radius of curvature 


shortening ratio = 2u,/(2aP/AE) 


The ratio 2u,/P in the elastic range was obtained 
from the slope of a straight line faired through the plot 
of measured values of shortening against load. The 
shortening ratio was computed from this by dividing 
by (2a,AE). The results are given in Table 6. The 
shortening ratio would be | if the reinforcement coun- 
teracted the effect of the hole completely. It would be 
5 for a circular hole without reinforcement in an infinite 
sheet. 

The shortening of the diameter of the hole at the mid- 
thickness of the sheet is plotted as a function of the load 
in Fig. 4a for panels 6, 14, and 8 and in Fig. 4b for panels 
The panels in each group differ only in 
It can be seen from 


1,.9,. and 7. 
the reinforcement of the hole. 
Fig. 4 and the values in Table 6 that, in both the flat 
sheet and the curved sheet, the diameter of the hole 
shortens the same amount until a load of about 40,000 
Ibs. is reached, which is near the maximum load for the 
flat panels. The reinforcements used on these panels 
were ineffective in reducing the shortening under load 
of the diameter of the hole. 

In Fig. 5, the shortening of the diameter of the rein- 
forcement is plotted as a function of load. The solid 
line is taken from Fig. 4 and is the curve of load versus 
shortening of the diameter of the hole at the mid- 
In both the curved and flat 


thickness of the sheet. 
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Shortening of diameter of reinforcement, in. 
(b) Panels with 10in. radius of curvature 
Figure 5.- Relation between shortening of diameter of reinforce- 
ment and load for flat panels and panels with 10 
inch radius of curvature 


panels, it is evident that the reinforcement fails to de- 
form with the sheet. The asymmetrical doubler plate 
shortens less than the symmetrical doubler plates, prob- 
ably due to the fact that the mid-thickness of the asym- 
metrical reinforcement is farther from the line of the 
load than the mid-thickness of the symmetrical ones. 


(2.3) Buckling and Bending of Panels 


There was no clearly defined critical or buckling load 
observed or measured during the tests of the panels. 
There was usually one half-sine wave buckle running 
the full length and width of the panel. 

The data showed the presence of bending strains in 
the unreinforced panels at loads near the initial load. 
The reason for the presence of this bending could not 
be determined. The bending strains were taken as 
one-half the difference between strain measured by 
gages placed back to back. At low loads, there was a 
linear relationship between the bending strain and the 
load; this relationship held throughout the elastic 
range for those panels with the small radii of curvature, 
and it held for a decreasing portion of the elastic range 
as the radius of curvature increased. The bending 
Strains for the unreinforced panels were expressed as the 
ratio of measured bending strain, ¢,,, at symmetrical 
positions on the panel to average elastic strain, «. = 
P/AE, far from the hole. If no bending was present, 
this ratio would be zero. 
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TABLE 7 
Bending Strain Ratios for Unreinforced Panels 
€y,t > oy 
-——Gage-—— 
Coordin- 
ates* ——Radius of Curvature (In.)——— 
x y 10 20 40 80 co 


1.75 0.00 0.45 0.44 0.41 0.41 0.38 
2.50 0.00 0.32 0.35 0.29 0.36 0.27 


5.50 0.00 O08 =... 0.10 0.20 0.20 
9.00 0.00 —0.02 —0.02 0:06 0.08 


0.00 3.50 0.06 0.05 0.04 —0.05 —0.07 
0.00 5.00 0.17 0.09 —0.15 —0.21 —0.08 


*x = distance from center of hole along centerline transverse 
to direction of loading, in. y = distance from center of hole 
along centerline parallel to direction of loading, in. 

t «,, = bending strain parallel to longitudinal centerline. 


The bending strain ratio is given in Table 7. It is 
plotted in Fig. 6 against distance from the center for 
gages on the transverse and longitudinal centerlines. 
From Fig. 6, it can be seen that the bending strain ratio 
approached a value of 0.5, or about one-sixth the me- 
dian fiber strain ratio, at the hole on the transverse 
centerline. The value of the ratio dropped rapidly as 
the distance from the center increased for the panels 
with small radii of curvature (panels 1-3) and more 
slowly for the panels with the large radii of curvature 
(panels 4 and 6). The data for panel 5 were not in- 
cluded, since they showed large differences in bending 
strains on the two sides of the hole. 

The scatter of values of the bending strain ratio on 
the longitudinal centerline 5 in. from the center is prob- 
ably due to variations in load distribution over the ends. 

The bending strains were not computed for the rein- 
forced panels, since the reinforcement and sheet did not 
act as a unit, and therefore values of ‘“‘bending strain”’ 
obtained from the gage readings on the sheet and rein- 
forcement opposite each other would have been mean- 
ingless. The plots of Fig. 7 for panel 12, in which the 
load is plotted against strain for the gages 0.10 in. from 
the hole on the reinforcement and sheet, indicate that 
the sheet and reinforcement did not act as a unit. 
The axial strain in the sheet is large compared to that 
in the reinforcement up to a load of 60,000 Ibs. Above 
60,000 Ibs., the reinforcement begins taking part of 
the load. 

The lateral deflections near the hole on panels with 
large radii of curvature were larger for the unreinforced 
panels than for the reinforced ones, indicating that 
bending of the panel had been reduced by reinforcing 
the hole. Differences in deflections at points far from 
the hole may be ascribed to experimental scatter. The 
deflections of panels 1, 7, and 9 with a 10-in. radius of 
curvature and with no reinforcement, symmetrical 
reinforcement, and asymmetrical reinforcement, re- 
spectively, are plotted in Fig. 8. It can be seen that 
the asymmetrical doubler plate is more effective than 
the symmetrical one in reducing the deflection near 
the hole. This is not surprising, since the asymmet- 
rical doubler plate has a larger moment of inertia 
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Fic. 6. Bending strain ratios in the direction of loading along the transverse and longitudinal centerlines for panels with no re 
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Figure 3 .- Effect of reinforcements on deflections of panels with a 
0-inch radivs of curvature; load, 144,6001b., maximum toads 
from 160,000 Ib. to 196,000 Ib. 


about the mid-thickness of the sheet and thus more 
readily resists bending. The panels with small radii 
of curvature showed only small deflections (deflections 
<0.5 thickness) in the neighborhood of the hole prior 
to failure, whereas the nearly flat panels showed large 
deflections (deflections >2 thickness) extending over a 
larger portion of the panel. For panels with 120-in. 
radii of curvature, the reinforcement reduced the de- 
flections at the edge of the hole by about 20 per cent 
from those obtained on the unreinforced panel; the re- 
duction was still appreciable along a line 1.5 in. from 
the outer edge of the hole. 


(2.4) Failure of Panels 


All of the panels were tested to failure except panel 
6. The load at failure and corresponding average 
stresses, based on the cross-sectional area far from the 
hole, are given in Table 1. 

The buckle pattern of one half-sine wave in the longi- 
tudinal direction and one half-sine wave in the trans- 
verse direction, except as modified in the neighborhood 
of the hole by the reinforcements, was characteristic 
of all the panels just prior to failure. The panels with 
radii of curvature of 80 in. to infinity retained this same 
buckle pattern at maximum load; the buckle pattern 
was obtained from deflection readings taken at maxi- 
mum load or just after maximum load had been reached. 
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The panels with the smaller radii of curvature, 10 to 40 
in., Showed a uniform pattern of a half-sine wave in both 
the transverse and longitudinal directions at loads prior 
to maximum load but showed an irregular buckle pat- 
tern at or just after maximum load. Either this ir- 
regular pattern was symmetrical on both sides of the 
longitudinal centerline or the deflections on one side of 
the longitudinal centerline were a mirror image of those 
on the other side. It is this difference in buckle shape 
that probably accounts for the maximum stresses being 
so much higher for panel 10 than for panel 2 and being 
lower for panel 11 than for panel 3. Panels 1 and 9 
had similar deflections before and after maximum load. 


(3) ANALYSIS 


(3.1) General 

The test data were analyzed in terms of the plane 
stress theory for the stresses and displacements in the 
neighborhood of a reinforced hole in an infinite flat 
sheet. The effect of curvature was neglected in the 
analysis, since no theory for curved sheet with a rein- 
forced circular hole was available. 

The following equations are given in reference 5 for 
the strain in a radial direction and the strain in a cir- 
cumferential direction, respectively : 


G 
= : 
E 


r? . a‘ a? 
—A(i+v) —6»B 7 —3C(1 + v)— — 2D — |cos 26 
a y4 r? 


){ra- ”) + K(1 +») —+ 
" 


(1) 
Ss f : . a? 
& = -)< Fil — v) — K(1 +7) =+ 
E r? 
: r? " a‘ a? 
A(1 +0) + 6B = + 3C(1 + ») — + 2Dv = |cos 26 
a? r’ r? 
(2) 
The following equation for the shortening of the radius 
of the hole in the direction of load is given by setting 


r = aand @ = 0 in the equation for the radial displace- 
ment [Eq. (14), reference 5]: 


Ua = (*=*\| ra —y— K(1 + y= 


A(1l+yv) —- 2B+ C(1 +p) + 2D | (3) 


The constants A, B, C, D, F, K can be determined from 
Fig. 3, reference 5. 


(3.2) Plane Stress Solution for Plate with Unreinforced 
Hole 
For a hole with no reinforcement, the ratio 7'/t of 
thickness of panel in reinforced region to thickness of 
panel in unreinforced region is 1, and the ratio 6/a of 
the outer radius of reinforcement to radius of hole is 1. 
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Figure 9 .- Median fiber 


for unreinforced 


The following values for the constants in Eqs. (1) to 
(3) are obtained from Fig. 3, reference 5: 
With 7/t = land }/a = 1, 


A=-' D=1 
B=0 F =1/, (4) 


C=-', K=-—'/, 


Substituting these values in Eqs. (1) and (2) gives 
the strain ratio in the direction of load along the longi- 
tudinal centerline, 6 = 0, 


e-E €; a er 
= = = 1 — 2.65 — + 1.95 (5) 
0.,/6=0 ae eer r* hi 


=) 


and the strain concentration ratio in the direction of the 
load along the transverse centerline (9 = 90°), 


fegke - _ a? a‘ ' 
( = = 1+ 0.05 — + 1.95 (6) 
o., 70=90° € @=90 eo ” 


The strain ratios in the direction of the load along the 
longitudinal and transverse centerlines were computed 
from Eqs. (5) and (6), respectively, with a = 1.488 in. 
They are plotted in Fig. 9 for comparison with the ex- 
perimental values for panels 1 to 6 given in Table 3. 
The agreement between experiment and theory is good. 

The shortening across the hole in the direction of load 
was computed for panels 1 to 6 as the ratio u,E/o,a 
from Eq. (3) and the constants in Eq. (4), with the 
result 


9 


Uk /o.a = 3 


panels. 


The experimental value of this ratio (Table 6) varies 
from 2.79 to 3.55. There is good agreement between 
the strain ratio at the hole and the shortening ratio, 
except for panel 5. 


(3.3) Plane Stress Solution for Plate with Reinforced 
Hole 


For a hole with the reinforcement shown in Fig. 1, 
the average value of 7'/t, the ratio of thickness of rein- 
forcement and sheet to thickness of sheet far from the 
hole, was 3.043, and the average value of b/a, the ratio 
between the outer and inner radius of the reinforce- 
ment, was 1.408. “1” and “2” 
indicate the unreinforced and the reinforced portions of 
the sheet, respectively, the following values for the 
constants in Eqs. (1) to (3) are obtained from Fig. 3, 


Using the subscripts 


reference 5: 


With 7 /t = 3.043 and b/a = 1.408, 
A, = —0.500 A, = —0.442 
B, = O B, = 0.040 
C, = 0.504 C, = —0.361 (7 
D= 0.421 D= 0.762 
Fi, = 0.500 F, = 0.302 
K, = —0.090 K, = —0.302 


If, in addition, Poisson's ratio v is taken as 0.», the 
strain ratio in the direction of load along the longitu- 
dinal centerline, 9 = 0, for the reinforced region (1.52 < 
r < 2.14) from Eq. (1), is 
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results for panels with a reinforcement on one side of 
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the symmetrically reinforced panels was 2.38 and 2.26, 
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and that for the panels with reinforcement on one side 
of the hole only varied from 2.46 to 2.86 (see Table 6). 
The reinforcement reduced the strain concentration at 
the hole but not so much as the theory indicated it 
should. The symmetrical doubler plate was more 
effective in reducing the strain concentration than the 
asymmetrical doubler plate. 

The elongation across the hole in the direction of the 
load was computed as the ratio u,//o,a from Eq. (3), 
the constants for the reinforced ‘‘2”’ region in Eqs. (7), 
and taking v = 0.3, with the result 


tat/o.a = 2.21 


The experimental value of this ratio, given in Table 
6, varied from 2.60 to 3.18 for those panels with asym- 
metrical reinforcement and from 2.78 to 2.92 for those 
with a symmetrical doubler plate reinforcement. The 
reinforcement was not so effective as the theory indi- 
cates in preventing the shortening of the hole under 
load, although the shortening was reduced somewhat. 
The shortening was less for the symmetrical doubler 
plate reinforced hole than for the asymmetrically rein- 
forced one. 

The reinforcement for panels 9 to 14 shortened less 
than the sheet, as can be seen from the values in Table 
6 and from Fig. 5. The reinforcements did not act asa 
unit with the sheet and, not being in line with the load, 
were not shortened so much. 


(4) CONCLUSIONS 


The tests indicate that, in the elastic region, the 
plane stress theory for flat sheet can be used to deter- 
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mine the median strain concentration and the shorten- 
ing of the diameter of the hole in a cylindrically curved 
panel with an unreinforced circular hole. The theory 
will also reliably predict strains in the sheet outside the 
reinforcement in panels with doubler plate reinforced 
holes. If the reinforcements are attached with only 
one row of rivets, the strains in the reinforcement will] 
be lower, and the strains in the sheet in the reinforced 
region higher, than those computed from the theory. 
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The Lift Distribution on Conical and Non- 
conical Flow Regions of Thin Finite Wings 


In a Supersonic Stream 


THEODORE R. GOODMAN? 
Cornell Aeronautical Laboratory, Tne. 


ABSTRACT 


A method is presented for finding the part of the pressure distri- 
bution which results in lift on wings in a supersonic stream. The 
linear theory isemployed. The method will solve for flow regions 
that are conical, as well as those that are not conical. The lead- 
ing edges may be supersonic or subsonic. When the trailing edges 
are subsonic, the Kutta condition is shown to be automatically 


satisfied in one particular case. 


INTRODUCTION 


— rWO BASIC LIFT PROBLEMS for a finite wing are: 
Given the lift distribution over a given plan form, 
find the twist and camber that will support such a load- 
ing. The other so-called direct problem is: Given the 
twist and camber on a given plan form, find the loading 
that the wing will support. 

The first of these two problems, the inverse problem, 
can be solved by direct integration. The second is 
more difficult, and only certain shapes have thus far 
yielded their solutions to the aerodynamicist ' 

Wings that have regions of conical flow may readily 
be solved by two possible techniques. One technique 
uses transformations in the complex plane; the other 
uses singular integral equations. These techniques 
may be used whether the leading edge is supersonic or 
subsonic. These methods are illustrated by Stewart,’ 
Brown,* and Heaslet, Lomax, and Jones. 

Evvard® has given a method that may be used easily 
when the so-called external flow regions do not interact, 
When they do 
For the 


whether the flow is conical or not. 
interact, the method becomes cumbersome. 
case of two subsonic edges that meet at a point, the 
external flow regions become infinitely interacting and 
the method fails completely. 

A method suggested by P. A. Lagerstrom (that of 
“canceling’’ the excess lift) is not under the restriction 
of conical flow. Nor does it have the restriction of 
noninteracting external regions. The method is, there- 
fore, powerful. This method is illustrated by Cole- 
man.°® 

The present paper offers a method of solving for the 
lift distribution in regions of a finite wing which is also 

Received December 15, 1948. 
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fairly general. The idea of canceling excess lift (or the 
principle of superposition, which is the same thing) is 
used, except that the excess lift is not canceled by 
superimposing elementary conical lifting surfaces. In- 
stead a method of supersonic doublets is used to carry 
out the cancellation. The technique is introduced by 
illustrative examples, starting with a simple problem 
The flow 


may or may not be conical, and the plan forms may have 


and ending with more complicated ones. 
subsonic or supersonic leading and trailing edges. 


THE IDEA 


The equation for the velocity potential in steady 
three-dimensional flow under the assumption of small 
perturbations is well known. It is 


(MI? — 1) (0°¢/0X*) = (0%y/O¥?) + (0%¢9/0Z*) (1) 


where 
M = the free-stream Mach Number 
X, Y,Z = the cartesian coordinates, where X is in 
the direction of the free stream 
re) = the perturbation velocity potential 


Suppose the problem were to find the twist and cam- 
ber of the wing for a given loading. Heaslet and Lo- 
max! have given the solution to this problem in terms 
of supersonic doublets. If u, v, w represent the Y, Y, 
Z perturbation velocities, respectively, and if the load- 


ing AC, is given for a wing lying in the Y, Y plane, then 


u  B?| f f Z ACA(X1, V1) dXid¥i 
U- 4n\|J J, [(X — Xi)? — BY — ¥)? - B27] 


(2) 
where 
B =VM-1 
AC, = Ap/q lower surface — Ap/g upper surface 
U = free-stream velocity 
T = the region on the wing subtended by the fore 
Mach cone of the point X, Y, Z 
gq = dynamic pressure [= (1/2)pU?] 
p = free-stream density 


The symbol means the ‘‘finite part of’’ (see ref- 
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erences 8,1). Now 
w/U = (0/0Z) f*.. (u/U)dX (3) 
The boundary condition is 


a(X, Y) = Lim (w/UV) (4) 
zZ—0 
where a is the local slope in Y = constant planes meas- 
ured with respect to the free-stream velocity vector. 

Thus, for a given AC,, the twist and camber of the 
wing may easily be determined. 

Eqs. (2), (3), and (4), however, may be thought of 
as an integral equation for finding AC, when a(X, Y) is 
given as pointed out in reference 1. It will be shown 
how this point of view may be used in finding AC, for a 
great variety of polygonal plan-form shapes. 

THE METHOD 

Before presenting the method, it is convenient to 
make a transformation; let 


A = Xy = XxX 

BY=y BY,1=y 

BZ = 2 BZ, = 2) 
Then Eqs. (2), (8), and (4) may be written 


O re . . ; — 
Lim 5 / dx / / K AC, (x1, v1) dxy dy, = 
sa 0 Z sai r 


a(x,y) (0) 


where 
. = B2/4e[(x — x1)? — (y — ym)? — 27] 
Since, from Eq. (5), the limit is taken as z > 0, 7 in 


the new coordinates will be taken as the region sub- 
tended by the fore Mach cone of the point x, y, 0. 


Tip of Rectangular Wing 

Consider a simple case, the tip of a rectangular wing. 
This problem may be easily solved by other methods, 
but it will serve to illustrate the present method, which 
may be used in more difficult cases. 


TICAaAtL 
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The problem is to find the pressure coefficient in the 
region between the tip Mach cone and the right edge of 
the wing (see Fig. 1). Suppose, temporarily, that the 
tip were nonexistent and that the wing were infinite 
Call the pressure coefficient under these 
The region + would then be re- 
1, and Eq. (5) 


spanwise. 
circumstances AC,,. 
gions (1) + (2) + (3) as shown in Fig. 
would be 


Lim 2 fa wf f K AC, 
¢ D+ 2+C6 


Assume that this problem has been solved and that 
other method) AC,, is 


(XY), V1) dx), dy, 


a(x, y) 6) 


from Eq. (6) (or by any 


known. 

Now return to the problem of finding the pressure 
coefficient in the tip region, between the Mach line 
from the origin and the x-axis. Call the pressure co- 
efficient in this region AC,,. Then the pressure co- 
efficient in region (1) will be AC,, by virtue of the fact 
that region (1) is not influenced by the tip. Eq. (5) 
for finding the pressure coefficient in region (2) is, there- 


fe re, 


Lim 2 fas ff K AC,,(%1, ¥1) dx, dy, + 
Lim 2 a SS KAC,( (x1, V1) stl = 


a(x, y) (7) 


The first term in this equation is identical with the 
left-hand side of Eq. (6) except that the region (1) re- 
places the regions (1) + (2) + (3). If, therefore, the 
first term in Eq. (7) is extended to include regions (2) 
and (3), then from Eq. (6) this will cancel a(x, y) 
which appears on the right-hand side of Eq. (7). This 
can be done provided AC,, is analytically continuous 
across the Mach line from the origin and across the x- 
axis. 

Eq. (7) now becomes 


> - —— ow 
Lim / dx / / K AC,,( 
¢—+ 0: OS — © \ )+@2)+(6) 

/ / K AC,,(x1, V1) dx dy, — 
| / / K AC,,(%1, Yi) dx dyy a a(x, vy) (8) 
2)+(3 f 


Now cancel the first term on the left with a(x, y). 
Also let f = AC,, — AC,,. Then Eq. (8) becomes 


iam of a 2) Af (x1, V1) dx, dy, = i 
im as dx | ff fC 90 a y1) dx, dy, (9) 


X1, ¥1) dx, dy, 5 i 


an 


wh 
rig! 
it i 
(on 


tair 


f(z, 


the 
of 
the 
‘ite 


ese 


ire 
ine 
CO- 
CO- 
act 


(9) 


(9) 


Lit? 


Differentiating with respect to z, taking the limit as 
z— 0, and then differentiating with respect to x, 


SS: K'f(x1, 1) dx, dy, = 
S Sey K' MC, (X1, 91) dx dy, (10) 
where 


K’ = 1/4x[(w¥ — 1)? — (y — y)?7]” 

It is seen that the lift that would occur to the right 
of the tip if the effect of the tip were not taken into 
account has been “‘canceled’’ and the effect of this 
canceling on the wing itself has been resolved into an 
integral equation for finding f. 

Now transform the integrals to the so-called Mach 


coordinates. Let 


MENT My = (Xi + ¥1)/2 
VW =m — My uy = (X%1 — V1) /2 
x =v+u v = (x + y)/2 
y=uv-u w= (x — ¥)/2 
The area element under this transformation is 
dx, dy, = (du, dv,)/2, and Eq. (10) becomes 
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Sf f(t, v1)duy do; 

| >) (u — mm) “*(v — 2) 
Sf AC,, (1, 11) duty dry 
| 3) (ua — 


Now, examine Fig. | again. 


ie (11) 
u;) *(v — v%) 
The equations of all 
the important lines are given in Fig. 2 in the u, v co- 
ordinate system. 

Putting in the limits, Eq. (11) becomes 





i dv, ( {"f(m, viduy ™ AC,,(t1, v1)diay 
3 f a fo 
0 w~— a)” mn (&%@ — %)° — 1 (u — u) f 


The solution to integral equations of this type is de- 
rived in the Appendix. Employing the notation of 
Whittaker and Watson,’ if the integral equation is given 


by 


1+ 


fix) = |f0* u(é)dt/(x — &) (0< p< 1) 


The solution according to the Appendix is 


—p sin wr / = f(x)dx 
l _ 
7 e @=-z) * 


Considering the quantity in brackets to be the unknown 


u\z) = 


and applying this solution to Eq. (12), there is obtained 


| | " f(u1, vr)duy / AC, (1, 01)duy 13) 
y= (1. 
| Ju (u%@ — m4)” —n (u — %)” 


where the ‘‘finite part of’ sign has been dropped on the 
right-hand side of this equation since it is finite. This, 
it is seen, amounts to dropping the outer integration 
(onv in this case). 

Applying the solution again to Eq. (13), there is ob- 
tained 


fle 1 J . du / * AC,,(t1, 01) duy 
4%) ==> os 1/o ~ SJ. 
2a Jn (2 —u)*SJ—-n (u — uy)” 


(14) 


(12) 


Interchanging the order of integration, this becomes 


;* 
f(z,%1) = — Qn J. AC,,(t1, v1) du, X 


* 


/ ’ du 
mn (u — m4) “(sc — u)” 


Using Peirce’s integral tables,® the first integration is 


(15) 


easily performed : 


f V2-—7, F kk AC, (ua, v1) dy 
jis, Vy *= — > - 
Tv —n % — Gis — M3) 
(16) 


Until now, no statement has been made concerning 
the twist and camber conditions under which AC,, was 
derived. Suppose the wing to be a flat plate at angle 


of attack, a. Then AC,, is a constant and equal to the 


two-dimensional Ackeret value—i.e., AC, = 4a/p. 
Then, 
) 2 4a , | 20; ir 
f(z,%1) = — tan | (174) 
7 wT OE \ a= & 


Let ¢ = uw, and transform back into x; and ¥; coordi- 
nates. Then, 


f(x1, v1) = —(da/Br) cos~! [—1 —2(y,/x1)] (18) 
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Therefore . 


pe = (4a/B) — Ua r83)\a _ cos~'[1 + 2(91/%1)]} 
(19) 


| =o ; = 
if, duy ry f (m1, v1)dv, 

: : 28 = 

| J (b/2)[(m + 1)/m] (u — uy)*”? m—b(v — 1)" 


As before, the outer integrations are identical. Heaice, 
the integrals on 7; may be equated. Eq. (21) becomes 


ee Saar ag = ~/ 
f(m, v)dv; i. AC,,(u1, 01) dv 
u,—b (v — v1)" R [(1 — m)/(1 + m)]ua w= 01)” : 


(22) 


_ 


| 





and the solution is 


l ‘ dy 
ee eo 
AC,, (1, 01) d01 


I. —b 
(23) 
[Cl — m)/(1 + m)]ur (v _— v1)" ° 


Interchanging the order of integration and integrat- 
ing, one finds, for AC,,, 
ssashdieaniesindsia 

fs : Vze-—-m+50 
AC,,(t1, 2) = AC, —- - Raden pd 


yA 


| 
{Ql — m)/(1 4 


Transform the variable of integration by letting 
t = yi/x1 = (v1 — m)/(%1 + 1) 
Also, let s = v; and transform back into cartesian co- 


T 
AC,, (u1,01)dvy on 
4a 2 = (24) 
m) ux (z oi 1)V uy, — b — vy 


Swept Forward Wing with Subsonic Trailing Edge 


Lett = y:/xi. Then 
AC,, = (4a/7B) cos~' (1 + 2#) 20) 


This is a well-known result. 


Delta Wing with Cutoff Tips 


To illustrate the method in another case, consider 
the delta wing with cutoff tips (see Fig. 3). 

It is proposed to find the loading in the region influ- 
enced by the cutoff tip where the pressure coefficient 
in region (1) (AC,,) is known. Since this is one un- 
known region superimposed on a known region, it is 
similar to the preceding problem and, with a similar 
definition for f, Eq. (11) holds true for this problem. 
Putting the limits in Eq. (11), the integral equation to 


be solved is: 


F — . 
I. duy | AC,, (m1, V1)d0 
yA ao amas 7? 
(b/2)[(m +1)/m] (U — 1) AS (Ll —m)/l+m]u (v — 2)” 


(21) 
ordinates. Then, 
j - Vy + i 
AC, (1, 1) = AC, (x1, 1) — ———— X 
T 
1/[1 — (x1 — y1)/6) 
(my — Va) f x 
— 98 
AC,,dt 
(25) 


V1 = ty, — xt)V —t(x1 — nm) — 001 — D 


For a flat plate wing at an angle of attack a (see 
references 2, 3, or 4): 


4am 1 


—— — =< (26) 
BE(m!') V1 om (y1/mx)? 


AC, (x1, ¥1) a 


or 


AC,,(t) _ a. _ os : 
BE(m") V1 — (t/m)? 

where E (m’) is the complete elliptic integral of the 
second kind and m’ = V1 —m*, Eq. (25), together 
with Eq. (26), is identical to the solution derived by 
Coleman in reference 6. 


The method will now be applied to solve a problem the solution to which has not been previously obtained. 
Consider the problem of finding the solution behind the trailing-edge Mach cone of a sweptforward wing (see Fig. 4). 


(20) 


ider 


iflu- 
ient 
un- 
t is 
ilar 
lem, 


n to 


(21) 


the 


ther 
| by 


ned. 
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This wing has subsonic trailing edges. It is generally accepted that the Kutta condition must be satisfied 
along a subsonic trailing edge—i.e., as one approaches the trailing edge AC, — 0. The Kutta condition 
for the present problem will be temporarily ignored. 

Eq. (11) may also be used for this problem. Substituting in the limits, the integral equation to be solved for this 


problem is 


ve dv, = f(m, m)dm 
<A i = 
c/2 (U — Vy) °A[Q + m)/(1 — m)] 2 — mc/(1 — m) (U — 1) 
| ; ro ey Es ar oe 
f dv, Tow on See Sy AC,, (ui, v1)duy - 
4 : (24) 
c/2 (9 — Bj)" Ja (% — %) °° 


Following a procedure similar to those in the preceding problems, there is obtained 


“ 


_ [ a> a <] 

xy + Vy y ~ Lar + on — me 

AC, (x1, 91) = AC, (%1, M1) — — Vm(e — 1) —- xX 
ve - . ‘ P : 0 


a x 


AC,,dt = 
—=—_— = (28) 
(ni = tx;)V1 + tV (m + t) (4:1 + yi) — mc(l + 2) 


It can readily be shown by direct evaluation of AC,, at the trailing edge that this solution satisfies the Kutta 
condition. Fora flat plate wing at angle of attack, a (see reference 4) : 


AC, (¥1, 1) = [2a BE(m')]¥ 2(1/x) (1 — V1 — m?)/[m — (91/x1)] 
or 


AC, (t) = [2a BE(m')|V¥ 2411 — V1 — m*)/(m — 8 


Wing with Small Aspect Ratio this equation, consider the derivation of Eq. (11). It 
was necessary to have a term that would cancel a so 
that the integral equation could be set down in such a 
way that the outer integration could be ignored. The 
same must be true here. It is necessary to find the 
term that will cancel a. For this purpose consider the 
wing shown in Fig. 5, assuming it were required to find 


Now consider a more difficult problem, one for which 
Eq. (11) doesnotapply. (See Fig. 5.) 

It is proposed to find the pressure coefficient in re- 
gion (5) knowing the pressure coefficient in the regions 
ahead. For this problem Eq. (11) no longer is valid, 
and a new equation must be derived. Before finding i piles oh ' 

the pressure coefficient in region (4). 


Then Eq. (5) becomes (sce Fig. 6): 
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" Now return to Fig. 5. The region 7 for finding the | 
"YY; pressure coefficient in region (5) differs from r for finding 
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-. (| a” a 
— a ae ~“ . K AC,, dx; dy. + 
——_—_— l = 
/ [ KAC,, dx; dy, + JS KAC,, dx; dy, 
ty 2 3 
IS KAC,, dx dy, =a (30) 


transform to Mach coordinates (where dx; dy, = 
du, dv;/2). Then, indicating the limits according to 


Fig. 6, 


Li 2 fe 1 ff xac l ie + 
Am — xs | . v 
0202 *, \| : p, GN, ¢ 1 
} oa : . f nu ma — 
Sf KAC,, du, dv, + fe duy J KAC,, dv, + 
3 |J 11 ul 
if du J Kae, do, =a (31) 
|/ IL I ; f 


the pressure coefficient in region (4) in only one respect. 
Consider Fig. 6. As one moves out along the line u, = 
u from the point (u, v), one crosses the line I before 
crossing the line III. From Fig. 5 it is seen that this 
is reversed in finding the pressure coefficient for region 
(5). 

It is now required to write Eq. (5) to find the pres- 
sure coefficient in region (5). But since it will be neces- 
sary to cancel a, first take the corresponding terms 
in Eq. (51) and then “‘correct’’ these so that the true 
region 7 for finding the pressure coefficient in region (5 
is left. Writing the respective terms in Eq. (31) and 
indicating the limits from Fig. 5, 


| a x ~~ 
Lim / dx / / K AC,, du, dv, + 
——taer-~ | 


- “ =; 
iff K AC,, du, dv, + (/ duy | K AC,, dv, + 
: 1 II] 
| u Sil "a fu ‘ my 
[ du; f K AC», ao.) + f du; 4 KAC,, dv 
JV I J 1 I 


— 


The terms in the above expresson will equal a since 
they correspond to the terms in Eq. (31). However, 
they do not represent the correct expression for finding 
the pressure coefficient in region (5). Eq. (5) for find- 
ing the pressure coefficient in region (5) is: 


ja f° cc 
Lim / dx J / / KAC,, du, dv, + 
——oson 4 , , 

lf KAC,, du, dv, + / / K AC,, du dv, — 

/ / i AC,, duy dv, + | | K AC,, diy day 4 
. 7 (6 4)-+ (5) +(6 

J | Kf du, dv, — / f K AC,, du; day 

K AC,, du, diy = a (0) 


where f = AC,, — AC,, and the assumptions of analytt- 
cal continuation are assumed to hold. 

The first five integrals in this equation are identical 
with those in expression (32), and, since they equal a 
from Eq. (31), Eq. (33) becomes 


o 4 ¢ | 
rs f _ dx J f Kf du, dv, — | 
J J ke du, dost = 


where g = AC,, — AC,,. Taking the derivative with 
respect to z, letting z — 0 and differentiating with 


Lim 
z— 0 


| 
' 


respect to x, this becomes 
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{ [ f duy dv, 
JSG (u — 4) (uv — 2%)” 


uv 


g duy dv, sae 
. SJ— * 3 (30) 
| Dt — f° w= hh) 


This equation is fairly general. From it one may 
find the pressure coefficient behind a Mach wave that 
is reflected from an edge provided the previous regions 
have been determined. For example, one could find the 
pressure coefficient in region (5) for the wing shown in 
Fig. 7 provided the solution were known in regions (:3) 
and (4). This problem has subsonic trailing edges, and 
it would be necessary to show that the solution satisfies 
the Kutta condition. 

As an example of the use of Eq. (35) the problem 
shown in Fig. 5 will now be solved. 

Placing the limits into Eq. (35) from Fig. 8 


| ' dy / 

b (uM — M4)? Ju—b (Vv — %) 
| du | ie 
| Jb (u — m)’* Jo 


and proceeding as before, there is obtained for a flat 


F(a, v1)dv% 


b 

g(%1,%1) dv 
(36) 
(v — 2) 


plate at angle of attack, a, 


oe 
cos 1+ 


ACp, ./ *=b/(xy— 4 — b) 
-¥v M1 +. b ia = V1) x 


AC,, AC,, + ) 


2(26-—m+y ] 
(XY) op V1) Z 


‘ — 1 
cos~! (1 + 2t)dt 
— (37) 


V1 —tin -m) V — te, — 9) —- 81-8 


where, for a flat plate at angle of attack, a: 


DISTRIBUTION ON WINGS 


IN A SUPERSONIC STREAM 371 


(AC,,/m) cos~'[(871 — m1)/(#%1 + u1)] — AC,, 


4 = 

AC,, = (AC,,/r)im — cos—' — [1 + (2y9:/x1)] — 
cos~! — [1 + 2(yv1 + b)/x]} 

AC,, sa 4a/p 


After some reduction and an integration by parts, 
this solution can be shown to be identical to that ob- 
tained by Coleman. ® 

Consider another problem (see Fig. 9). 

It is proposed to find the pressure coetlicient in re- 
The procedure will be similar to that used 
First 


gion (7). 
in finding the pressure coefficient in region (5). 
formulate Eq. (5) for finding the pressure coefficient 


in region (5). The regions are shown in Fig. 10. 
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ae ng si “e 
Lim : ae dx 1 : K AC,, dx; dy, + 
/ fa K AC,, dx dy, a j [ K AC, dx dy, a+ 
- « arma ——— — > 
| / K AC,, dx; dy, + / / K AC,, dx; dyay =a 


(38) 





Transforming to Mach coordinates and indicating 
the limits according to Fig. 10, 


ss eS 2 ee 

Lim dx - K AC,, dv, du, + 
s—-0202 J-—~o@ | 1 

[fo svi ye wT Tr aa 
f dv / K AC,, du, + / / K AC,, dv, du, + 
VI J Iv JG 

lS: Fil : ay ar a ~~ 
I a, f K AC, du, + | duy I K AC,, divi ¢ = a 
| Jvi VII I Vv f 


(39) 


Now, write Eq. (5) to find the pressure coefficient in 


region (7). From Fig. 9 the terms corresponding to 


those in Eq. (39) are 


; 1 a *y j | . 7 ; . 
Lim dx ) K AC,, dv, du, + 
PA >0 202. oe . 1 
*1 evil . 
( / dv J K AC,, du, + 
JVI J Iv 
{0 VII — oe 
/ dv, / AC. du, ) + / / K AC,,dv,du, +4 
mp JIV . a Gs 
- TI =a S , 
/ dv, J K AC,, duy / duy / BAC. dv ( 
JVI J vil Jt JV 


(40) 
Eq. (5) for finding the pressure coefficient in region 


(7) is (see Fig. 9): 
Lim 4. / dx / / KAC,, doy du, 
2-0 2 Oz — 2 ( - 7 il 
[/ K AC,, dv, du, — | / / K AC,, dv, du, + 
Ii. K AC,, dv, du, + | If K AC,, dv, du, + 
3 le 4)-+(6)+ (8 
ff K AC,, dv, du, + yi KAC,, dv, du, — 
5)+G m 


SS KAC,, dv, du, + ee Kf dv du, =a 
{ 6+ | : f 


(41) 
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where f = AC,, — AC,, and the usual assumptions of 


‘ 
analytical continuation are assumed to hold. 

The first six integrals in this equation are identical 
with those in expression (40), and, since they equal a 
from Eq. (39), Eq. (41) becomes 


m7 ie oo ee 
Lim dx 4 | K AC,, dv, du, — 
Oz J — « | J JG i 


2-0 


// K AC,, dv, diy + | If Kf dv du! 0 


i) 


Differentiating with respect to 2, setting s equal to zero, 
and differentiating with respect to x as before, this be- 


comes 


Fi ) 
f dv, duty 


JJ 4) 
. 7) (% — %&) 
TS. 4+@ (u — mm)" (v — % 
/ / AG. dv du 13 
J JS@(u — mm)" (v — 9)” 


3) is general for a region behind two Mach 
Thus, for the wing of 


A (Uv — Be) 


AC a dv, du, 


Eq. (4 
waves that cross on the wing. 
Fig. lla (the same wing as in Fig. 7), for example, the 
pressure coefficient may be found in region (7), provided 
it is known in regions (2) and (4). 

Likewise, referring to Fig. 11b, the pressure coefii- 
cient may be found in region (7), provided it is known 
in regions (2) and (4), by using Eq. (43). But Eq. 
(43) may also be used to find the pressure coefficient in 


re: 


th 
fre 


th 
su 


sic 


LIFT DISTRIBUTION 


region VII of Fig. 11b provided it is known in regions 
II and IV. Also, since Eq. (35) may be used to find 
the pressure coefficient behind a Mach line reflected 
from an edge, it is seen that by successive application 
of Eqs. (85) and (43) the pressure coefficient in all the 
regions of the wing shown in Fig. 11b (and lla) may 
Only regions (2), (3), (4) require different 
handling. But for regions (2) and (3), Eq. (11) is 
valid, and for region (4) another simple integral equa- 
tion may be derived in a similar manner. (It is, of 
course, necessary to check the Kutta condition on the 


be found. 


wing shown in Fig. 11a.) 


THE METHOD IN GENERAL 


|) Write Eq. (5) for a point (x, y) or (4, v) as 
though it were for a point in a region where the pres- 
sure coefficient is known. 

(2) Write Eq. (5) for the point (x, y) or (#, v) con- 
sidering it to be in the unknown region. 

3) Subtract the equation obtained in step (1) 
from the equation obtained in step (2), canceling a. 
An integral equation for f remains where f is the differ- 
ence between the pressure coefficient in the adjacent 
region (which is known) and the pressure coefficient in 
the unknown region. 

1) Differentiate with respect to zs, set s = 0, and 
differentiate with respect to x. 

(5) If the outer integration on all remaining inte- 
grals is the same (in Mach coordinates), equate the inner 
integrations. 

(6) Apply the solution given in the Appendix to the 
resulting integral equation. (If the solution given in 
the Appendix does not apply, one may have to resort 
to numerical methods. ) 

It is to be noted that the term that cancels a (step 1) 
must be derived from a region immediately preceding 
the unknown region or else the outer integration on the 
resulting integral equation will not be the same. 

The method given in the present paper for solving for 
the pressure distribution in new regions always requires 
a knowledge of the first region. This offers no par- 
ticular difficulty, however, since the first region can 
usually be solved by the methods given in references 
2,3, 4,5 

Th 
gions on all possible plan-form shapes. 
to define what may be called Mach convex and Mach 


e present method will not solve readily for all re- 
It is of interest 
concave polygons. It will be recalled from plane geom- 
etry that a convex polygon is a polygon for which no 
straight line exists which intersects more than two of its 
Also 
from plane geometry it will be recalled that the in- 


sides. Concave polygons are all other polygons. 
ternal sides of a concave polygon are the sides that form 
an internal angle of the polygon which is greater than 
180°. 
po'ygon for which a Mach line exists which intersects 
more than two of its sides and which has two adjacent 


A Mach concave polygon will be defined to be a 
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Fig. 12 


internal subsonic sides, both of which lie within the aft 
Mach cone (or both within the fore Mach cone) of the 
point which forms their apex. Mach convex polygons 
are all other polygons. 

Now, provided the pressure distribution in all the 
‘first’ regions on the wing can be found by some other 
method (and there may be more than one ‘“‘first’’ re- 
gion as evinced by the two tips of a sweptforward wing), 
then the method presented in this paper will readily 
solve for the pressure coefficient in all the other regions 
on the wing in terms of a quadrature, provided the plan 
form is of the Mach convex type. 

Mach concave polygons will always offer some diffi- 
culty. The pressure coefficients in two unknown re- 
gions may interact—i.e., a second unknown region may 
appear within the fore Mach cone of a point in the first 
unknown region and vice versa, thus giving rise to two 
simultaneous integral equations for the pressure coefli- 
cients in both regions. The solution given in the 
Appendix will not apply, and the integral equations 
may have to be solved numerically. This difficulty 
appears in the solution to the symmetric wing with a 
subsonic trailing edge when the two branches of the 
trailing edge come to a point. 

Another difficulty that may occur in solving Mach 
concave plan forms is that the Mach lines that define 
each new region may reflect in such a manner as to 
give rise to an infinite number of regions before the plan 
form is completely covered. [Mach convex plan forms 
may also have an infinite number of regions (such as 
Fig. lla) but never before the plan form is completely 
covered. Thus, for these plan forms, the total lift can 
always be found to within as close an approximation as 
required.] This phenomenon: occurs in the swept- 
forward wing with a subsonic leading edge when the 
two branches of the leading edge come to a point. 
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It is always possible by a slight modification to make a 
Mach concave plan form into a Mach convex plan form 
and thus avoid the mathematical difficulties inherent in 
Mach concave plan forms. If the two adjacent subsonic 
edges that are causing the difficulty are not permitted to 
come to a juncture but a small supersonic edge ts allowed 
to replace the juncture, the plan form becomes Mach con- 
vex, and then the pressure coefficient in all the regions 
can be solved in terms of a quadrature. 

This may be seen with the aid of Fig. 12. 
shows the sweptback wing with a subsonic trailing edge. 
A small supersonic edge has been introduced at the 
juncture of the two branches of the trailing edge (shown 
as a dotted line), and the regions induced by the junc- 


This figure 


ture are also shown as dotted Mach lines. 


APPENDIX 
Suppose the given integral equation is 


F(x) / *  u(&)dé 
: Tr 
a (f — £) 


(0< n< 1) 


Integrate the right-hand side of this equation by parts. 


Then, 


u(&) 


‘a u'(é)dé 
u(x — &)* Ja (x — €)* 


f(x) = = 


Taking the finite part, 


f u(a) f u’(&) dé 
(x) = — - 
u(x — a)” bJa (x — €)? 


where it is assumed that u (a) exists. 
This integral equation for u(£) may be written 


f u'(t) dé 
Ja (x — &)* (Xx 


f u(a) 
— uf(x) — 
— @)" 

This is an integral equation of the Abel type. Its 
solution is given in Whittaker and Watson.’ The solu- 
tion to the integral equation 


. * u'(&) dé 
f(x) = : — (O< w< 1) 
we 1 = S&S) 


is given by Whittaker and Watson to be 


sin war : f(x) dx 
T Ja lo=— x)! . 


u(a) + 


u(s) = 
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Applying this to the problem under consideration, the 
solution is seen to be 


: ; psinur [2 f(x) dx 
u(z) = ula) — ——— 
» u 
T a (g — xX) f 


sin ur : ‘ dx ' 
u(a) ; 
T 4 (x — a)*(g — x) 74 


The transformation x = (s — a)y + a changes the 


second integral into 


. ad | 
sin ur, dy 
u(a) —*- 
T o y(1 — y) * 


sin ur 
aa ula) (w)T(1 — pw) = (a 


from a well-known identity. Therefore, the solution 


for u(z) is 


u(z) = —(u sin wr /7) So’ f(x) dx/(z — x) 
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A Method of Least Square Adjustment for 
Measured Wing Deflections 





DALE R. BRANCHFLOWER* 
Northrop Aircraft, Inc. 


SUMMARY 


he dynamic characteristics of an airplane wing are predicted 
from calculations based on the theoretical bending-torsion modes 
The validity of these predictions depends upon 
For this pur- 


of the wing 
experimental verification of the theoretical modes. 
pose measured wing deflections are obtained from ground vibra- 
tion tests 

rhis paper applies the principles of least square adjustment to 
the measured modes in such a manner that a best fit is obtained 
Furthermore, the measured deflections 
Probably the 
most important feature of the method given is that the coefficients 


for the kinetic energy. 
ire expressed as a series of theoretical deflections. 


in the series are a numerical measure of the similarity of theoreti- 


cal and experimental modes 


INTRODUCTION 


fpr COUPLED, BENDING-TORSION MODES of an air- 
plane wing calculated from design data are nor- 
mally compared with measured modes obtained by 
ground vibration tests. The degree of similarity of the 
two sets of modes indicates the degree of validity of flut- 
ter characteristics and dynamic landing loads calculated 
from theoretical modes. In ground vibration tests, 
deflection amplitudes and phases are measured at a 
number of points while the wing is harmonically ex- 
cited to resonance. The bending and torsional dis- 
placements so determined do not ordinarily form smooth 
plane curves when plotted against a spanwise coordi- 
nate. The conventional method of fairing smooth 
curves through the plotted points is subject to several 
objections: For instance, faired curves cannot be du- 
plicated independently; further, there is no assurance of 
an optimum fit, nor is there any measure of the dis- 
sinilarity of the faired and the computed curves. The 
first two objections may be met by a suitable applica- 
tion of the ‘“‘Method of Least Square Adjustment’’' to 
the measured data, and the third objection is obviated 
by introducing the theoretical modes as expansion func- 
tions, as shown in the following discussion. 

WING THE Mass DISTRIBUTION 


COORDINATES AND 


An arbitrary, spanwise, reference axis of the wing is 
assumed to form a plane curve having a continuous 
first derivative. Any point on the wing is located by 
the distances s and 7 measured along the reference axis 
and normal to it, respectively. The coordinate s is 

Received July 6, 1948. Revised and received January 7, 1949. 
\erodynamicist Engineer. 





measured outboard from the root chord, and r is posi- 
tive for points aft of the reference axis. 

The wing area is partitioned into quadrilateral areas 
by drawing straight lines normal to the reference axis 
at the points s; and intersecting them with spanwise 
lines, denoted by the index 7. Thus, the lattice points 
of the network correspond to the numbers 7;;(s;). 

The mass distribution of the wing is approximated by 
a set of concentrated masses m;;(s;) at the points 7;;(s;). 


LEAST SQUARE ADJUSTMENT OF MEASURED WING 
DEFLECTIONS 


The vertical deflection in the kth theoretical mode is 
Since the functions f;;“ form a 


represented by f;;. 
, may be 


complete set, the sth measured mode, Z;, 
represented by the series 


Zy™ = Var fy™ (1) 
k 


Let w,;“? be a set of measured deflections for the lat- 
tice points. If the true deflections are given by Eq. (1), 


define weighted residuals by 
(s) | ge ae 7 (8 9 
ej = Vmi (wy — Zu) (2) 
The ‘“‘most probable’ deflection function is that for 
which the sum of the squares of the residuals is a mini- 


mum. 
the coefficients a,“. The “‘normal equations’ follow 


This minimum is obtained by a proper choice of 


from the condition 


(0/da,) ) S>(e,)?2 = 0 (3) 
(aj \ 


Substituting from Eqs. (1) and (2) into Eq. (3) and 
simplifying, 


(s) a P l % ae r (kK 
da, muffs ( = Vom iywy Fi (4) 
l (i,j ) iJ 


Let h; and g;“ represent the bending of the refer- 
ence axis and the twisting around the reference axis, 
respectively, in the kth theoretical mode. Then, 

fy = h® + rag (5) 

Therefore, 


muff iy fy = 
1.3 


fm hh, D+ Slag + hg) + 


Jig} (6) 
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where 


mn; = Yomi; 
J 

S; = Somiry (7) 
J 


‘e 9 
myn 
j } 


The theoretical coupled modes satisfy the orthogo- 


Nw 
II 


nality condition? 
DY {mh On +. Sih Og + hig, ] ES 
1 
Figg} = K 6x1 (8) 


where 
K, = Dom{(hi)? + 2S Ag; + Silgi)? 
; (9) 
. _ flk=l 
"VOR +1 


Substituting from Eqs. (5), (6), (8), and (9) into Eq. 
(4), 


a, = (1/K,y) om ywy fy 
J 
= (1 K,) omyw (hk, + ryg;) (10) 


tJ 


Numerical values of the a,‘’, computed from Eq. 
(10), are substituted into Eq. (1) to obtain the adjusted 


deflections. 


EQUATIONS FOR THE BENDING DEFLECTION, TORSIONAL 
DEFLECTION, NODE LINE LOCATION, AND KINETIC 
ENERGY 


From Eqs. (1) and (5), the bending and_ tor- 
sional deflections in the sth measured mode are, respec- 


tively, 
k 


a? = Dare | 


4, = eet 


(11) 


Let 7; be the normal distance to the reference axis, 
corresponding to the point s; on the axis. Similarly, 
let Z,;° represent the vertical deflection of the wing at 
r; Both r; and Z;° 
point s; The points for which Z,“ = 0 lie on the node 


are continuous variables for each 
lines of the wing. Therefore, 
Z:9 = Yah + 74%q,) = 0 
p 


and 


ri) = +e (12) 


If w, is the measured mode frequency, the kinetic 


energy, /,, is given by 


EB = 2) w,2 om (24; )? (13) 


t 
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From Eqs. (1), (6), and (9), 
Tmy(Zy)? = Ly Daa, Ofy! 
i i) 1% \ 
= Yaa) Ymuyfyf yj | 
kl \ 


44,7 
= (a, ’ )*K, 
k 
Therefore, 


E. = ({1/2)e » * (a;,“° )*K, 14) 
k 


The coefficients 2E,/w,? are required in both flutter 
analyses and dynamic load calculations based on meas- 
ured modes. Eq. (14) presents a simple way of obtain- 
ing the numerical values of these coefficients. 


DISCUSSION 


Replacement of the distributed wing mass by a set 
of point masses introduces an error in the section mo- 
ments of inertia, as given in Eq. (7). This error be- 
comes negligible if the number of point masses is suf- 
ficiently large. 

The introduction of the weighting function, V m,,, 
into the residuals, Eq. (2), is responsible for the occur- 
rence of orthogonal sums in the left member of Eq. (4). 
The physical significance of the weighting function in 
the least square adjustment process is to make a best 
fit for the kinetic energy, expressed in Eq. (14). 

A useful consequence of the orthogonality condition 
is as follows. The calculation of each coefficient a," is 
independent of the calculation of the remaining coef- 
Therefore, higher approximations to the de- 
, may be found without duplicat- 


ficients. 
flection function, Z;;“° 
ing previous work. Additional terms are simply added 
to the sum in Eq. (1). 
Note that, if w,;° = fi; 
every coefficient vanishes except that for which / = k, 
Therefore, for any other 


’ is substituted into Eq. (10), 


and in this case a,“ = 1. 
deflections, w;;“, the coefficients a,“, 1 + k, are meas 
ures of the dissimilarity of the actual deflections and the 
theoretical mode f;;“. 

Theoretical coupled bending-torsion modes are usu- 
ally calculated with respect to an “elastic axis.”’ The 
elastic axis is defined as that axis for which elastic cou 
pling terms vanish from the differential equations of 
If the elastic 


axis is chosen as the reference axis in the least square 


motion, leaving only inertia coupling. 


adjustment procedure, comparison with theoretical 
modes is facilitated. 

The analysis in the preceding sections has assumed no 
control surface deflections. If the measured displace- 
ments include rotation of control surfaces, some assump- 
tion must be made regarding the relative displacements 
with respect to the wing chord. If the assumed control 
surface displacements are subtracted from the measured 


(Continued on page 383) 
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Reinforced Monocoque Cylinders 


John E. Duberg 
Aeronautical Research Scientist, Structures Research Division, 
Langley Aeronautical Laboratory, N.A.C.A., Langley Field, Va. 


May, 1949 


— READER IS in substantial agreement with the material 
presented in the paper ‘Shear Stress Concentration and 
Moment Reduction Factors for Reinforced Monocoque Cylinders 
Subjected to Concentrated Radial Loads” by N. J. Hoff, Vito 
L. Salerno, and Bruno A. Boley (JOURNAL OF THE AERONAUTICAL 
ScieNcES, Vol. 16, No. 5, May, 1949). I feel, however, that the 
paper may be misinterpreted, as it now stands, regarding the 
action of the structure as a whole. 


When an isolated ring in a structure is loaded, there is a reduc 
tion in the maximum bending moment in the loaded ring below 
that computed in the usual way and, correspondingly, a large 
increase beyond that predicted by conventional analysis in the 
maximum shear stress in the adjacent sheet. There is also a 
transfer of bending moment from the loaded ring to the unloaded 
rings. It can be shown that for a structure with many rings, at 
a given azimuth, the sum of the bending moments in the un- 
loaded rings approaches the difference between the elementary 
computed bending moment and the actual bending moment for 
For a cylinder with infinitely many rings, the 
From a practical point 


the loaded ring. 
sum is precisely equal to the difference. 
of view, only those rings in the immediate vicinity of the loaded 
ring participate in the bending action. The number of rings in 
volved depends on the relative stiffness of the rings. A glance 
at the charts given in N.A.C.A. T.N. No. 1310 (reference 8 of 
the subject paper) indicates that for stiff rings almost all the dif- 
ference in bending moment is shared equally by the two rings 
idjacent to the loaded ring. 
larger values of A and 4/B), more rings participate in the action 


As the rings become less stiff 


and this becomes less true 


This discussion suggests that some care should be used in com 
puting bending moments in rings of monocoques when more than 
one ring is loaded and, in particular, when many rings are simi 
larly loaded. A particular example is one of the monocoque 
cylinders (No. 3) tested and stress analyzed in N.A.C.A. War 
time Report L-66 (reference 5 of the subject paper). The ele 
mentary solution gives a maximum bending-moment coefficient 
in any loaded ring equal to 0.239Pr. When each of the four 
rings is loaded separately, the maximum bending moments in 


the loaded rings are: 


Ring 1 2 3 i 
Moment 0.207Pr 0.167 Pr 0.160Pr 0.147Pr 
Decrease 0.032 0.072 0.079 0.092 


When all the rings are similarly loaded, the maximum bending 


moments in the rings are: 


Ring 1 2 3 4 
Moment 0.248Pr 0.233Pr 0.212Pr 0.171Pr 
Decrease (—0.009) Pr 0.006Pr 0.027Pr 0.068 


It can be observed that for this particular case it is even possible 
for a maximum bending moment to exceed the value given by 
elementary considerations. Ring one, for which this is true, is a 
ring at the tip and is supported from one side only. 

Similarly, some attention should be directed to the shear 
stresses when more than one ring is loaded. The primary or 
elementary shear stress arises from the relative tangential dis- 
placement of points on the loaded ring (temporarily considered 
rigid and therefore circular) with respect to the adjacent rings. 
This distribution is sinusoidal, and the division of it between the 
bays to either side depends on the supporting conditions of the 
cylinder as a whole. The major portion of the secondary effect 
on the shear stress distribution that arises from ring bending is 
caused by the relative tangential displacement of points on the 
loaded ring and points in the same azimuth position in the un- 
loaded rings. When many rings are loaded similarly, there is no 
large relative tangential displacement, due to ring bending, be- 
tween the loaded rings, and therefore the shear stresses are sub- 
stantially those given by elementary theory. If conditions of 
symmetry eliminate the warping of the cylinder cross sections, 
that portion of the secondary shear strain arising from this action 
is also eliminated, and the shear flows between rings will be pre- 
cisely that given by ordinary engineering theory. Under such 
conditions the secondary effect will be, for the most part, confined 
to the region between loaded and unloaded rings. 

A summary of the shear stress analysis for the same cylinder 
of reference 5 gives some indication of these effects. From ele- 
mentary theory the maximum shear flow in a circular cylinder 
loaded by a radial force is 0.318(P/r) and occurs 90° from the 
load. When each ring of the test cylinder was loaded separately 
the maximum shear flow occurred in the adjacent inboard bay 
approximately 45° away from the load and had the following 


values: 
Bay 1 2 3 iY 
Shear Flow 0.434(P/r) 0.483(P/r) 0.532(P/r) 0.675(P/r) 


0.116 0.165 0.214 0.357 


Increase 
When all rings are loaded the elementary solution yields the fol- 
lowing set of values for shear flow: 


Bay 1 2 3 


4 
Shear Flow 0.318(P/r) 0.637(P/r) 0.955(P/r) 1.273(P/r) 


The maximum values obtained by considering the secondary ef- 


fects are: 


Bay 1 2 3 1 
Shear Flow 0.338(P/r) 0.646(P/r) 1.012(P/r) 1.540(P/r) 
Increase 0.020 0.009 0.057 0.267 


It is clear that for this loading condition on this cylinder there 
has been a considerable reduction in the secondary effect 
except in the bay immediately adjacent to the rigid root connec- 
tion. 

I believe the foregoing indicates the danger of misapplication 
of the factors derived for the case of only one loaded ring to 
problems involving several loaded rings. 


Authors’ Reply 

N. J. Hoff, V. L. Salerno, and Bruno A. Boley 
Polytechnic Institute of Brooklyn 

May, 1949 


— AUTHORS WERE very much interested in reading Dr. 
Duberg’s valuable remarks on the effects of loading simul- 
taneously several rings in a reinforced monocoque cylinder 
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This problem was not treated in the original paper. There is 


complete agreement between Dr. Duberg and the authors re- 
garding the theoretical aspects of the subject. On the other 
hand, calculations carried out at the Polytechnic Institute of 
Brooklyn have shown that the deviations from engineering 
theory remain extremely large even in the case of several similarly 


loaded rings, provided the parameter A has a high value. The 
following numerical example, based on a three-field cantilever 
having A = 10’? and B = 500, will illustrate the point in ques- 
tion: 


1) First ring alone loaded 


Ring or Ringfield (Bay) l 2 3 
. 

Maximum moment 0.081 Pr Negligible Negligible 
Conventional maximum moment 0.239Pr 0 0 
Decrease 0.158Pr be 
Maximum shear flow 2.91P/r 2.12F 1.89P/r 
Conventional maximum shear 

flow 0.318P/r 0.318P/r 0.318P/r 
Increase 2.59P/1 1.80P/r 1.57P/r 
Il 1/1 the ec rip ¢ 
Ring or Ringfield (Bay) l 2 3 
Maximum moment 0.096Pr 0.089 Py 0.082Pr 
Conventional maximum moment 0.239Pr 0.239Pr 0.239P; 
Decrease 0.143 Pr 0.150P; 0.152Pr 
Maximum shear flow 2.18P/r 4.06P/r 6.18P/1 
Conventional maximum shear 

flow 0.32P/r 0.64P/r 0.96P/r 
Decrease 1.86P/r 3.42P/r 5. 22P/s 
Ill) Second ring alone loaded 
Ring or Ringfield (Bay) 1 2 3 
Maximum moment Negligible 0.079Pr Negligible 
Conventional maximum moment 0 0.239Py 0 
Decrease its 0.160Pr 
Maximum shear flow Negligible 3.01P/1 2.23P/r 
Conventional maximum shear 

flow Negligible 0.32P/r 0.32P/r 
Increase Negligible 2.69P/r 1.91P/r 
1V) Third ring alone loaded 
Ring or Ringfield (Bay) 1 2 3 
Maximum moment Negligible Negligible 0.073Pr 
Conventional maximum moment 0 0 0.239Pr 
Decrease 0.166Pr 
Maximum shear flow Negligible Negligible 3.46P/r 
Conventional maximum shear 

flow 0 0 0.32P/r 
Beewemme ll —“‘“‘“Cs™s™sSCSC‘*SSCC ee 8 ene 3.14P/r 


It appears therefore that the quantitative conclusions reached 
by Dr. Duberg are valid only for small values of A—that is, for 
unusually heavily reinforced fuselages of small diameter. In 
fact, the cylinder of the NACA series referred to had A = 1815 
and B = 15.9. Since in modern transport planes the value of A 
is usually between 107 and 10%, and that of B between 250 and 
400, the numerical example given here is more applicable to ac- 
tual design than the N.A.C.A. test cylinder 
ample shows only small changes in the values of the stress con- 


The present ex- 


centration and moment reduction factors because of the simul- 
taneous application of loads to several rings 


Dr. Duberg’s physical interpretation of the behavior of the 
cylinder and his observation that in an infinitely long cylinder 
the sum of the bending moments in all the rings at any given azi- 
muth is equal to the conventional bending moment constitute 
The last statement is ap- 
parently approximately correct even for a short cylinder of only 
three bays when A = 1815. When the value of A is of the 
order of magnitude of 10’ corresponding to modern transport 
fuselages, the cylinder must have a great many rings between the 
loaded and rigidly fixed ends in order to make this quantitative 
conclusion hold even approximately. 


significant additions to the theory. 
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On Harmonic Motion at Supersonic Speeds 


John W. Miles 
Department of Engineering, University of California at Los Angeles 
February 23, 1948 


— LINEARIZED EQUATION for a time dependent velocity 
potential (@) in a supersonic, inviscid flow may be written 


l 2@ o\ 
V %4(x, y, 2, t) = l + (x, ¥, 
c? Ox ol 


where U is the free-stream velocity and c is the sonic velocity in 


the free stream. It is generally desired to solve this equation 
subject to the prescription of the velocity (i. e., Ag) at a surface 
that is approximately (in accord with the linearization) parallel to 
the flow 

If the harmonic time dependence exp(7zwt) is assumed, a con- 


venient transformation of Eq. (1) yields 


(x, y, 2,t) = (x, 9, exp|i(wt — Ax 9 
V y2*v = B2(O*p/0x?) + (w/Bc)*¥ 

N = [1 + (1/8?)] 7?(w/Be } 

B = [(U?/c?) — 1 5 


where V , a plane normal to the 


free stream 


2 is the Laplacian operator in 
Now if (w 


potential equation for steady flow; 


‘Bc)*? is neglected, Eq. (3) reduces to the 
accordingly, if the steady 
flow solution for the given boundary conditions is known, a first 
order (in powers of the frequency) approximation to the non 
steady flow problem may be constructed with the aid of Eq. (2 
Retaining only the first order term in A in the exponential, the 


result may be written 


g(x, y,2,t) = Fe dS(é, 7, €) [1 — tA(x — &)] &X 
S 


G(x —£,¥y— 9,2 — S)weé, 9, §, 0) + OA 6 


where (é, 7, ¢) is a point on the boundary surface, w is the pre 
scribed perturbation velocity at this point, and G is the Green's 
function for the problem at zero frequency and is, therefore, inde 
pendent of ». The pressure coefficient is given by Bernoulli's 
theorem as 


C,(x, y, 2, 2) = — p(x, y, 2, t)/(pU?/2) = 
2/0 “) 7 
a mm) x Vv ‘ 

Ll’ \ox [ 


As a simple example, the mid-chord moment on a two-dimen 
sional, flat plate airfoil, which is executing a pitching motion 
about its mid-chord, will be calculated. If x is measured from 
the leading edge, positive downstream, z is measured positive 
down, and the chord is two units of length, the displacement of 


the airfoil may be taken as 


z(x, t) = (x — 1) exp(iot 8 
and the down-wash is given by 
Oz Oz lw 7 
w(x,t) = + =/1+—_ (x — 1) } exp(tol 9 
ot Ox l 
The Green’s function for the steady flow problem is simply 
: 1 (1/8 (x > 0 
G(x) = l(x) =? 10 
f { 0 (x < O 


Substituting w and G in Eqs. (6) and (7) to obtain c,(x), the re 
sulting moment coefficient is 


1 > ; 
Cna * ; ,’ (x — 1)2c,(x)dx = 
2A\} 2(3 + 28?) jii- ‘ 
il - + 0(? lb 
38 6B? 1 + p* 


which, after allowing for the difference in notation, is in agree- 
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ment with the first term in the power series expansion of the 
flutter results.' 

The appropriate Green’s function for a quarter infinite wing, 
having the leading edge x = 0 and the port edge y = 0, is found 
to be 


_/ i l 1/2 (am —1/3 5 9 
G(x, y) = -» (x = ¥) l(y) l(x — y) (12) 
i re 

This result may then be used to obtain approximate pressure 
distributions on oscillating rectangular wings for which BA.R. is 
greater than two, and the lift and moment coefficients computed 
will be valid for BA.R. greater than one. If Cne is calculated 
for such a wing, the imaginary part agrees with the sum of the 
derivatives Cng and Cng given by Harmon.? [The writer has 
also obtained an exact solution to Eq. (1) for a rectangular wing. } 
In the case of a slender body of revolution, the potential is es- 
sentially independent of both Mach Number and frequency to 
the approximation justified by the original linearization and is 

given by the well-known Munk approximation,? viz. 


U' cos 0 


g(x, r, 0) = R2(x) w(x) (13) 


Yr 

where y and @ are polar coordinates measured from the x axis and 
the vertical, respectively; x is measured from the nose (assumed 
pointed); and R(x) is the radius of the cross section. The result 
is, of course, dependent on the frequency insofar as w is concerned 
e.g., Eq. (8)], and the pressure is further dependent on the fre- 
quency by virtue of the differentiation with respect to time, cf. 
Eq. (7 This latter point was overlooked by Beskin in a calcu 
lation of the damping moment on an oscillating body of revolu- 
tion.! 

It is felt that the foregoing results are of practical importance 
in flutter work, inasmuch as the most serious flutter problems— 
e.g., 1° torsional instability—appear to occur at low frequencies. 
Nevertheless, it must be remarked that the approximation here 
issumes small (w/c) and not simply small (w/c). 


Miles, J. W., The Acrodynamic Forces on Oscillating Airfoil at Supersonic 
Speeds, Journal of the Aeronautical Sciences, Vol. 14, No. 6, p. 356, Eq. 
A2), June, 1947 
?Harmon, S. M., Stability Derivatives for Thin Rectangular Wings at 
Supersonic Speeds, N.A.C.A. T.N. No. 1706, November, 1948. 
Munk, Max, The Aerodynamic Forces on Airship Hulls, N.A.C.A. T.R 
No. 184, 1923 
‘ Beskin, I Damping Forces and Moments on a Body of Revolution at 
Supersonic Velocities, Consolidated Vultee Aircraft Corporation, CVAC- 
DEVF Memo BB-5, May 31, 1946 


Gas Flow in Ducts 


Frank W/. Barry and C. Darby Fulton 

Research Assistants, Department of Mechanical Engineering, Massa- 
chusetts Institute of Technology, Cambridge, Mass. 

March 10, 1949 


_ IS A DISCUSSION of the paper ‘‘On the Addition of Heat 

to a Gas Flowing in a Pipe at Subsonic Speeds’ by Joseph 
V. Foa and George Rudinger (JoURNAL OF THE AERONAUTICAL 
Sciences, Vol. 16, No. 2, pp. 84-94, February, 1949). 

The subject of gas flow in ducts, or aerothermodynamics, has 
been developed during the last few years. Simple heating has 
received the greatest attention because it is basic to the ram-jet, 
gas turbine, and rocket. That even this special case is not easy 
is shown by the errors, both philosophical and mathematical, 
that have been made by competent writers. 

rhe philosophical error of Bailey in assuming that the de 
frease in total pressure implies irreversibility was discussed by 
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one of the writers previously. Samaras, in looking for a special 
manner of heating by which that decrease could be avoided, de- 
duced incorrectly from his Eq. (21) that the total pressure re- 
mains constant if heating is carried out at constant Mach Num- 
ber. That equation assumes an indeterminate value when the 
condition of constant M is inserted in it, a circumstance that ap- 
parently did not come to his attention. The correct relation for 
the ratio of total pressures as a function of total temperatures, for 
constant Mach Number, is (using the authors’ notation) 
P2/P, = (Ty PA 27M? 

The simple rule for minimizing the loss in total pressure is to add 
each increment of heat at the lowest possible Mach Number. 
More generally, the treatment of Hawthorne and Shapiro shows 
all the factors that influence the total pressure. 

The most frequently made error, and at the same time the most 
absurd, was the deduction that combustion stops when the Mach 
Number reaches unity or that heat cannot be added in this re- 
gion. It was never reasonable to think that gross equations of 
state, continuity, momentum, and energy could reveal chemical 
reaction rates. It was also unreasonable to deduce that, while 
the gas could not be heated according to a stationary observer, 
it could be heated according to a moving observer. The neces- 
sary degree of freedom permitting escape from this absurdity is 
found, of course, in the boundary conditions. 

The authors’ heavy emphasis upon boundary conditions is a 
good antidote. They have considered the following two cases: 
(a) constant m, 7), and also pp, when 1/2 < 1.0, and (b) constant 
P;, T;, and also p2 when AM, < 1.0; while other writers have 
considered, in considerably more detail, the case (c) constant 
M, nd 7. (Still other cases exist.) Unfortunately, the authors, 
upon obtaining some different results, imply that the other 
writers are generally in error. In so doing they fail to recognize 
what is valid and what is invalid in the previous papers, and they 
also reveal a narrower view of the whole problem than would 
have been expected of them. The authors have discussed 
transient phenomena resulting from the sudden addition of heat 
at a section in order to show how the inlet and exit conditions are 
altered by the heat transfer. They have also compared two 
flows, one being adiabatic, between the same boundary condi- 
tions and have shown how the inlet and exit states are a function 
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Fic. 1. General plot showing 44,, M2, AT/T;, and m curves 
for heating of air in a duct with constant inlet stagnation tem- 
perature and reservoir pressure. 
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Fic. 2. Showing drop in static temperature near end of 
duct as total heat addition is increased, when inlet conditions 
are held constant. Heat added uniformly alongduct. y = 
1.4, M, = 0.7. Back pressure must be lowered as heat is 


increased. 





x OR 


of the total amount of heat transfer and of the Mach Number of 
the initial adiabatic flow. Other writers have discussed only 
steady, established flows and have considered the change of state 
of a bit of gas as it traverses the tube. Most of them have 
chosen case (c) because it is analytically simple and because all 
other cases can be derived from it, and they have arrived at the 
appropriate boundary conditions as an end result. While errors 
have been made in failing to recognize the adjustment of the 
boundary conditions, the general method of considering the 
change of state of the gas as it traverses the tube is not only 
valid, it is indispensable. It is necessary to combine this ap- 
proach with a proper accounting for the boundary conditions in 
order to solve the whole problem. 

Fig. 1 is presented as an aid in obtaining a broader picture of 
what happens when heat is transferred to air flowing frictionlessly 
inapipe. For adiabatic flow the conditions would be represented 
somewhere along the diagonal line AT/7T; = 0. 
conditions when heat is flowing may be found at the intersection 
of the appropriate solid A7/7; line with the dashed line m = 
constant which is being followed. The ratio of total inlet pres- 
sure to exit static pressure may be found from the dotted lines. 
For case (b) the conditions when heat is flowing may be found 
at the intersection of the appropriate solid A7/7; line with the 
dotted line P;//) = constant which is being followed. The mass 
flow relative to the adiabatic mass flow may be found from the 
This plot shows that, as more heat is added, 


For case (a) the 


dot-dash lines. 
Mz increases for case (a) and decreases for case (b) and that M, 
decreases for both cases. The two sets of mass-flow curves are 
only applicable to the respective cases for which they were plot- 
ted. The other curves apply to any case, such as case (c) for ex- 
ample. This plot shows not only how the inlet and exit states 
change with the total heat addition but also how the state of the 
gas changes as it traverses the tube with a given total heat addi 
tion. The MM, = 1 line is an envelope of the P:/p9 curves; the 
value of P;/po at a point on the M2 = 1 line may be greater than 
or equal to the value of P;/po for the curve coming to that point 

While the authors have stressed the importance of boundary 
conditions, they have neglected to mention in their summary or 
conclusions the boundary conditions used. Their statement that 
the static temperature increases continuously even when > 
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1/+/y is completely incorrect for case (c), as many writers 
have shown. Fig. 2 shows how the static temperature varies as 
one increases the total heat addition, at the same time adjusting 
the back pressure or inlet pressure to maintain the inlet Mach 
Number constant. The authors attempted to show by the use of 
entropy that such a decrease of static temperature would be qa 
violation of the Second law of Thermodynamics. The proof 
could have been convincing only if a perpetual motion machine 
of the second kind had been shown to result. The authors have 
detracted from an otherwise excellent paper by misleading and 
incorrect statements in the summary, which is the part of the 
paper having the widest circulation. 

Four more references are included in this discussion to supple- 
ment the nearly complete bibliography of the paper: 

Edelman, G. M., and Shapiro, A. H., Tables for Numerical Solution of 
Problems in the Mechanics and Thermodynamics of Steady One-Dimensional! 
Gas Flow Without Discontinuities, Journal of Applied Mechanics, Vol. 14 
No 4, pp. 344-351, December, 1947. 

2? Shapiro, A. H., and Edelman, G. M., Tables 
Problems in Compressible Gas Flow with Energy Effects, Journal of Applied 
Mechanics, Vol. 15, No. 2, pp. 169-175, June, 1948. 

>’ Dusinberre, G. M., Hall, N. A., Rudnick, P., Shapiro, A. H., and Haw 
thorne, W. R., discussion and authors’ closure of reference 16, Journal of 
Applied Mechanics, Vol. 15, No. 2, pp. 179-181, June, 1948. 

4 Keenan, J. H., and Kaye, J., Gas Tables, Tables 30-53; 
Sons, New York, 1948 


for Numerical Solution of 


Tohn Wiley & 


Oscillating Airfoils at Mach Number 1 


N. Rott 

Institute for Aerodynamics, Swiss Federal Institute of Technology 
Zurich 

February 28, 1949 


- MAY BE OF INTEREST that for nonstationary flow finite solu- 
tions are obtained by linearized theory even for M = 1, in 
spite of the strong simplifications involved by linearizing in this 
case. This may give a further indication of the importance of 
nonstationary phenomena at sonic speeds." ? 

The basic solution at 14 = 1 for harmonic oscillations of a 
point source in three-dimensional space is given by*® 


yg = (A/x)expt —('/2)ik[x (y? + 2*)/x}} 1 


where the factor e*”* is omitted, k = v/a, and a = velocity of 
sound. For a source in a plane the following expression is read- 


ily obtained: 


¢ = (B/Vx)exp[—(#/2)ik(x + y?/x)] 2 
For an oscillating airfoil of infinite span in the plane y = 0, 


leading edge at x = 0, the y-component of the velocity on the 
surface is given in the form v(x)*’‘e. The potential is obtained 
by superposition from Eq. (2) and will fulfill boundary conditions 


at y = O putting 


~ 

* [ l dé : 

g=- o(x — E)exp| — = sk(E + y*/E) 7 3 
0 L ~ V 2riké 


The domains above and below the 


in the upper domain y > 0. 
airfoil are independent ahead of the trailing edge, thus by anti- 
symmetrical continuation of the potential [Eq. (3)] in the domain 
y < 0, the solution for the oscillating airfoil of zero thickness is 


obtained. (Symmetrical continuation would give the ‘‘pulsing’ 


plate.) 

Ackeret, J., and Rott, N., Strémung von Gasen durch ungestaffelte Pro- 
filvitter (The Flow of Gases Through Straight Cascades), Schweiz. Bauzet- 
tung, 1949, pp. 40-41 and 58 61 

2 Gardner, C. S., Ludloff, H. F., and Reiche, F., Drag of an Airfot ic- 


clerated Supersonic Flight, Journal of Applied Physics, p. 117%, 1945 ec 
also Aviation Week, p. 17, January 31, 1949 
2 Rott, N., Das Feld einer raschbewegten Schallquelle (The Field of 


ing Sound-Source), Mitt. Inst. f. Aerod. ETH, No. 9, 1945 
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To prove Eq. (3), the velocity v is calculated: 


oy . ie : : ; eae dt 
= = tky vo(x — t)exp}| — = h(E + 3°/2) | — 
Oy 0 “ Sv 2ntké 


or introducing 7» = y?/é, 
. sie i. = dn : 
= tk vo(x — y?/n)exp| — — tk(n + y?/n) a (4) 
SJ y?/x ~ V 2rikn 


Eq. (4) gives for y = 0 


. 


H < dn 
= tk v(x)exp | — — tkn ta = v(x) 
0 , 2 V 2rikn 


by the limiting value of the Fresnel integral ((!/2)kn = s): 
al os, 1 
] =-(1—*) = -;7 
Jo V2xs 2 V2 
The pressure 
pb = —p\(O¢/ot) + a d¢/dx)] 
is given on the upper surface of the airfoil by 
; f° > i 7 1 ° » vs 
p = pa [¢tkvo(x — £) + vor"(x — éE)Jexp] — = tke | X 
ZV ~ 
dé : E.. 
: + pavo(0)(2rikx)— /*exp | — = tkx (5) 
V 2rikt 2 
In the case of stationary flow—i.e., k = 0—p becomes infinite 


with k But for k ¥ 0, the pressure remains finite (except 
forx = 0), so that finite forces and moments are obtained for air- 
foils of zero thickness. 

Lift L and moment M (about mid-chord) have been calculated 
for a rigid airfoil (chord c) plunging (ye*’*) and pitching about 
the mid-chord (ae*”*). The result may be expressed in terms of 
the Fresnel integral of the reduced frequency w = (!/2)kc. Pro- 
visionally, the limiting values for small w should be given. Omit- 


ting again e*”* 


VW 4 l : 
, = yr —(2yo/c) - V-tw + ao X 
12 Vr ( 


1 7 
(. ze + 20 \ ww) (4) 


The factor 4/7 indicates a phase angle of 45°. The moment 
about the point x = ('/3)c (= 1/3 chord) vanishes as w > 0. 

If iv and accordingly 7k and iw are replaced by a real positive 
value (i.e., the Fresnel integral replaced by the error function), 
the same results hold for exponential divergence. 


On the Oscillating Rectangular Airfoil at 
Supersonic Speeds 


John W. Miles 


Department of Engineering, University of California at Los Angeles 
March 4, 1949 


i he PRESSURE DISTRIBUTION On a quarter infinite, zero thick- 

ness supersonic airfoil, having an arbitrarily prescribed dis- 
tribution of down-wash which exhibits the harmonic time de- 
pendence exp(iw!) can be found by using Fourier transforms to 
formulate the boundary value problem as a Wiener-Hopf type 
integral equation,! for which an exact solution may be given. 
' ' Miles, J. W., Transform and Variational Methods in Supersonic Acro- 
tynamics, Journal of the Aeronautical,Sciences, Readers’ Forum, Vol. 16 
No. 4, pp. 252-253. April, 1949 
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The result is valid for a rectangular wing where the Mach lines 
from the leading edge wing tips do not intersect on the wing. 

The expression for the pressure distribution due to an arbitrary 
distribution of down-wash over the wing is rather complicated, 
but the integration of the pressure over a spanwise strip of a rec- 
tangular airfoil effects a considerable simplification. Thus, if 
the down-wash is constant over the span, the lift coefficient for a 
spanwise strip x on a rectangular airfoil of leading edge x = 0, 
trailing edge x = 2, and aspect ratio A.R. is given by 


2 OL (x, t) 
orl, t) =|) = 4(M?-1 
p x 


g(x) = (2 a it) }[Jo(xx) — (2xA.R.,)~! sin (xx)Jexp X 
(—2ix Mx); (2) 
k = (wb/l (3) 
« = [kM/(M? — 1)} 1) 
AR, = (MW? — 1)? AR. (5) 


] re) . 
a(x,i) = — — + ik} s(x, t (6) 
U \ox 


where U is the free-stream velocity, ./ is the free-stream Mach 
Number, & is the usual reduced frequency parameter (based on 
the semichord 5), A.R., is the ‘‘effective aspect ratio,” and z is 
the displacement of the airfoil above its equilibrium position. 
This result differs from the two-dimensional solution only in the 
addition of the sin (xx) term, so that the three-dimensional cor- 
rection of the two-dimensional results? is a relatively simple 


matter. For the limiting case of steady flow (k = 0) and con- 
stant a, Eq. (1) reduces to 
Cp(x) = 4a (M? — 1)~'/? [1 — (2A.R.)~1x] (7) 


in agreement with the well-known Busemann result. 

Similar results for down-wash distributions with other span- 
wise weighting factors may be obtained, and it appears that the 
required integrals will generally be elementary, with the excep- 
tion of the integrals of Jo(xx), which have already been evalu- 
ated.? It is also of interest to remark that,.whereas the results 
(1) and (2) are valid only for A.R.g > 2, the various lift and 
moment coefficients obtained by integration are valid for A.R., 
ay. 

Numerical computations based on the above results (including 
flutter derivatives and gust response) are being undertaken in the 
Ballistics Division of the U.S. Naval Ordnance Test Station, 
Inyokern, Calif., and will be published, together with details of 
the analysis, in the near future. 

2 For a complete résumé of the two-dimensional results see Karp, J. N., 
Shu, S. S., Weil, H., Biot, M. A., Aerodynamics of the Oscillating Airfoil in 
Compressible Flow, F-TR-1167-ND, HQ, AMC, Wright Field, Dayton, 
Ohio, 1947; also F-TR-1195-ND, 1948. 

2 Busemann, A., Infinitesimale kegelige Uberschallstromung, Jahrbuch der 


Luftfahrtforschung, Vol. 7B, pp. 105-121, 1943. 


Formulas for the Theoretical Lift of Flat Wings 
at Supersonic Speeds 


Ellis Lapin 
Douglas Aircraft Company, Inc. 
March 3,.1949 


HILE THE THEORY of flat lifting wings at supersonic speeds 
has already been treated extensively, a summary of formu- 
las is convenient because of the somewhat lengthy and tedious 
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Summary OF FORMULAS FOR THE THEORETICAL LIFT OF FLAT WINGS : 
Supersonic LEADING EdaES 
at apis [Roa eho ae es sla ae Raa? PRR AAT 
te ae {Sr cos moe — Re mn cos" ine a (F sige! FOR m, =00 * 
ae ans (ne (ln mim ot - ae eerste +[Be4d - eta int St as (5-£&)} FOR m=m ‘3) 
qe am ("as len~ a +0, 008 nae +[ eee ~ Head] os sete « br pass yom + ara  } FOR m,=-m (4) 
Ba RE a [mo wv SS] a et ete ts cn Santee | 
= Safer {5 so se cos Banga) med ets ie: sine, FoR m=o0 (6) 
gt» Fs |i ae om (alo es) as ie oo eal 
Bae = aie | SLs Sas ma “el I a noe 2a gE FD coe mans 
7 ar Sm som Sussonic leap, ING EDGES — ” ' 

> om boda -or 3] . 

an aye" -co8'S} + a, fata} FOR m, = 00 (10) 

ae > ate pV ~ 2] + ay z} FOR m=m () 

it tg S (med, | ede a - ‘+ a(n] FOR m.=-m tia 

i =- ee g m)cos' + [ Seg on} WHERE: (3) 

“— Bebe} 8 rete | AND WHERE: 
488 jeomeet a sr ta| FoR m,=00 whERE m4’2)= red ee oa -a)- ad (4 
i= APPS, [cine & + (‘sae FOR m,=*tm WHERE (15) 
@, 
G'im oe [+ ats] FOR m=m | 

IN FORMULAS ABOVE WHERE £ OR F OCCURS, UPPER SIGN 1S TO BE USED WHEN M>+1 (SWEPT BACK TRAILING EDGES) 

AND THE LOWER SIGN IS TO BE USED WHEN m,<-1 (SWEPT FORWARD TRAILING EDGES). | 











integrations necessary for any particular wing. The formulas 


presented below were derived using the superposition of conical 
flows procedure as given by D. Cohen in N.A.C.A. T.N. No. 
“The Theoretical Lift of Flat Swept-back Wings at Super- 
e, with certain exceptions, extensions of 
of the 


1555, 
sonic Speeds,” and ar 
equations contained within that paper. 
notion that the lift of a flat wing may be obtained from the 
“basic” lift, which is the lift the finite wing would have if the 
side edges did not influence the situation, to which is added a tip 
Except for certain combinations of 
as noted below, 


Use is made 


induced lift correction. 
geometry and Mach Number, 
theoretically correct lift. 

The restrictions defining the class of wings to which these for- 
mulas apply are: The wings are flat, symmetrical about a plane 


this gives the 


normal to the wing, and trapezoidal; the side edges are parallel 
to the free-stream direction; the leading edges are swept back 
but may be subsonic or supersonic; the trailing edges may be 
swept back or swept forward. In the event the trailing edges are 
swept to an extent such that they are subsonic, an additional 
correction, described in T.N. No. 1555, must be made before the 
true lift can be obtained. Another restriction on the validity 
of the tabulated formulas is that, except for the cases m: = ©, 
the AL/ga formulas become invalid for those combinations of 
Mach Number and geometry such that the leading-edge tip 
Mach cones cross the wing centerline. For the cases m: an 
additional correction, not given here, is necessary when the 
leading-edge tip Mach cones cross the opposite tip. These limi- 
for most ‘‘practical’’ wings, the 


tations are not too serious, since, 





3) 


(4) 


(5) 


(9) 


(10) 





(i) 


(\4 


(15) 


— 


a 


allel 
back 


y be 





5 are 
onal 
> the 
idity 
is of 
tip 
an 
the 
imi- 
the 


READERS’ 


Mach Numbers at which they occur are generally close to unity, 
a region in which the applicability of linearized nonviscous super- 
sonic flow theory is dubious. 

The nomenclature is that of 7.N. No. 1555, where: 


ly) = basic lift (lift of entire wing ignoring the influence of 
the side edges) 


AL = tip induced lift correction 

’ = free-stream dynamic pressure 

a = wing angle of attack, rad. 

Wf = free-stream Mach Number 

re) = VM*?*-1 

x, ¥ = coordinates taken in chordwise and spanwise directions, 
positive aft and to the right, respectively, with origin 
at wing vertex 

ml = B(dy ‘dx)leading edge 


= B(dy dx) trailing edge 
= £6 (slope of ray from wing vertex to tip trailing edge) 
8 (slope of any ray from wing vertex) 
semispan 
E = complete elliptic integral of the second kind with modu- 
) 
lus (V 1 — m? 
The formulas may be made nondimensional through division 
by the wing area, S, and that part of the coefficient given by 


1s°/S may be replaced by 
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Examination of the lift formulas shows that B(dCy,/da) and 
8(AdC,/da) are functions of the parameters m;/m, a;/m, and 
m. There are in preparation charts of these functions in which a 
wide range of Mach Numbers and geometry are covered, by 
means of which the variation of lift curve slope with Mach Num- 
ber or geometry readily may be traced. 





Measured Wing Deflections 


(Continued from page 376) 


displacements, the present method becomes applicable 
to the derived data. 
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